#define _USE_MATH_DEFINES
#include <iostream>
#include <vector>
#include <string>
#include <cmath>

using namespace std;

const double a = 1;
const double b = 3;

/*
const double a = 1;
const double b = 3;
const double a1 = 1;
const double a2 = 1;
const double A = 2 * M_E;
const double b1 = 1;
const double b2 = 1;
const double B = 2 * M_E * M_E * M_E;
const double eps = 0.003;

string yx = "e^x";
double y(double x) {
    return exp(x);
}

string qx = "sin(x)";
double q(double x) {
    return sin(x);
}

string px = "cos(x)";
double p(double x) {
    return cos(x);
}

string fx = "exp(x) + cos(x) * exp(x) + sin(x) * exp(x)";
double f(double x) {
    return exp(x) + cos(x) * exp(x) + sin(x) * exp(x);
}
*/

const double a1 = 1;
const double a2 = 1;
const double A = 6;
const double b1 = 1;
const double b2 = 1;
const double B = 30;
const double eps = 0.1;

string yx = "2x^2";
double y(double x) {
    return 2 * x  * x;
}

string qx = "sin(x)";
double q(double x) {
    return sin(x);
}

string px = "cos(x)";
double p(double x) {
    return cos(x);
}

string fx = "12x + 6x^2cos(x) + 2x^3sin(x)";
double f(double x) {
    return 4 + cos(x) * 4 * x + sin(x) * 2 * x * x;
}

vector<double> TridiagMatrixAlg(const vector<vector<double> > & arrA, int n) {
    double h = (b - a) / n;

    vector<double> fi;
    for (int i = 0; i < n - 1; i++) {
        fi.push_back(f(a + i * h));
    }

    vector<double> m;
    vector<double> k;
    for (int i = 0; i < n - 1; i++) {
        double xi = a + i * h;
        m.push_back(-2 + h * p(xi));
        k.push_back(1 - h * p(xi) + h * h * q(xi));
    }

    vector<double> c;
    c.push_back((a2 - a1 * h) / (m[0] * (a2 - a1 * h) + k[0] * a2));
    vector<double> d;
    d.push_back(k[0] * A * h / (a2 - a1 * h) + fi[0] * h * h);

    for (int i = 1; i < n - 1; i++) {
        c.push_back(1 / (m[i] - k[i] * c[i - 1]));
        d.push_back(fi[i] * h * h - k[i] * c[i - 1] * d[i - 1]);
    }

    vector<double> y(n + 1);
    y[n] = (b2 * c[n - 2] * d[n - 2] + B * h) /
        (b2 * (1 + c[n - 2]) + b1 * h);
    for (int i = n - 1; i > 0; i--) {
        y[i] = c[i - 1] * (d[i - 1] - y[i + 1]);
    }
    y[0] = (a2 * y[1] - A * h) / (a2 - a1 * h);

    return y;
}

vector<double> fDiffMethod(int n) {
    double h = (b - a) / n;

    vector<vector<double> > arrA(n + 1);
    vector<double> fi(n + 1);
    for (int i = 0; i <= n; i++) {
        arrA[i].resize(n + 1);
        for (int j = 0; j <= n; j++)
            arrA[i][j] = 0;
    }
    arrA[0][0] = a1 * h - a2;
    arrA[0][1] = a2;
    fi[0] = A * h;

    arrA[n][n - 1] = b2;
    arrA[n][n] = b1 * h - b2;
    fi[n] = B * h;

    for (int i = 1; i < n; i++) {
        double xi = a + (i - 1) * h;
        arrA[i][i - 1] = 1 - p(xi) * h + q(xi) * h * h;
        arrA[i][i] = p(xi) * h - 2;
        arrA[i][i + 1] = 1;
    }

    for (int i = 1; i < n; i++)
        fi[i] = f(a + (i - 1) * h) * h * h;

    return TridiagMatrixAlg(arrA, n);
}

int main() {
    cout << "y'' + " << px << "y' + " << qx << "y = " << fx << endl;

    int n = 4;
    bool norm = true;
    vector<double> y2 = fDiffMethod(n);
    do {
        vector<double> y1 = y2;
        n *= 2;
        y2 = fDiffMethod(n);
        norm = true;
        for (int i = 0; i < y1.size(); i++) {
            if (fabs(y2[i * 2] - y1[i]) > eps)
                norm = false;
        }
    } while (!norm);

    cout << "y = " << yx << endl;
    cout << "n = " << n << " eps = " << eps << endl;
    for (int i = 0; i < y2.size(); i++) {
        cout << "y = " << y(a + i * (b - a) / n) << " resY = " << y2[i] << endl;
    }
}