#include <iostream>
#include <cmath>
#include <functional>
#include <vector>


using namespace std;

typedef std::vector<double> vec;
typedef std::vector<vec> Matrix;

Matrix Hein(double a, double h, double ya, double dya, double p, double q, std::function<double(double)> fx)
{
    Matrix Y(4);
    std::function<vec(double, vec)> f = [&](double xx, vec yy) -> vec{
            vec temp ={ yy[1], fx(xx) - p * yy[1] - q * yy[0]};
            return temp;
    };
    const int m = 10;
    vec curPt =  { ya, dya };
    h /= m;
    for (int i = 0; i <= 3 * m; ++i)
    {
        double x = a + i * h;
        if (i % m == 0) Y[i / m] = curPt;

        vec fi = f(x, curPt);
        vec newPt = { curPt[0] + h * fi[0], curPt[1] + h * fi[1] };
        vec fii = f(x + h, newPt);
        newPt = {
            curPt[0] + h / 2 * (fi[0] + fii[0]),
                curPt[1] + h / 2 * (fi[1] + fii[1]),
        };
        curPt = newPt;
    }
    return Y;
}

Matrix Solve(double a, double b, int n, double ya, double dya, double p, double q, std::function<double(double)> fx)
{
    double h = (b - a) / (n - 1);
    Matrix Y(n);
    vec F(n);
    std::function<double(double, double, double)> f = [&](double x, double y, double dy) -> double{
        return fx(x) - p * dy - q * y;
    };
    Matrix StartPts = Hein(a, h, ya, dya, p, q, fx);
    for (int i = 0; i < 4; ++i)
    {
        double x = a + i * h;
        Y[i] = StartPts[i];
        F[i] = f(x, Y[i][0], Y[i][1]);
    }

    for (int i = 4; i < n; ++i)
    {
        double x = a + i * h;
        Y[i] = 
        {
            Y[i - 1][0] + h * Y[i - 1][1] + h * h / 360 * (323 * F[i - 1] - 264 * F[i - 2] + 159 * F[i - 3] - 38 * F[i - 4]),
                Y[i - 1][1] + h / 24 * (55 * F[i - 1] - 59 * F[i - 2] + 37 * F[i - 3] - 9 * F[i - 4])
        };
        F[i] = f(x, Y[i][0], Y[i][1]);
    }
    return Y;
}


int main()
{
    double p = -3, q = 2, A = 4, B = -1, a = 0.15, b = 0.85, h = 0.05;
    int n = std::roundf((b - a) / h ) + 1;
    std::function<double(double)> fx = [&](double x) -> double
    {
        return x*x - 7*x;
    };
    Matrix temp = Solve(a, b, n, A, 0.13585390768417982, p, q, fx);
    for (int i = 0; i < n; ++i)
    {
        cout << "y(" << a+i*h << ")" << " = " << temp[i][0] << endl;
    }
}