#include <iostream>
#include <math.h>

#define EPS 0.001

using namespace std;

bool checkRunge(double* y1, double* y2, int s) {
    double max = 0;
    for (int i = 0; i < s; i++)
        if (fabs(y1[i] - y2[2 * i]) > max)
                max = fabs(y1[i] - y2[2 * i]);

    return (max < EPS);
}

double p(double x) { return -1. / 2; }

double q(double x) { return 3; }

double f(double x) { return 2 * pow(x, 2); }

double ai(double x, double h) { return 1 - h * p(x) / 2; }

double bi(double x, double h) { return -(pow(h, 2) * q(x) - 2); }

double ci(double x, double h) { return 1 + h * p(x) / 2; }

double di(double x, double h) { return pow(h, 2) * f(x); }

int main() {
    double a1 = 1, a2 = 2, A = 0.6, b1 = 1, b2 = 0, B = 1, a = 1, b = 1.3, h, x;
    int n = 3;
    int saverN = n;
    int i;
    double* y1, *y2, *pp, *qq;
    do {
        double b0, c0, d0;
        double an, bn, dn;
        h = (b - a) / n;
        b0 = -(h * a1 - a2);
        c0 = a2;
        d0 = h * A;
        an = -b2;
        bn = -(h * b1 + b2);
        dn = h * B;
        y1 = new double[n + 1];
        pp = new double[n + 1];
        qq = new double[n + 1];
        i = 0;
        pp[i] = c0 / b0;
        qq[i] = (-d0) / b0;
        for (x = a; i < n; i++, x += h)
        {
            pp[i + 1] = ci(x, h) / (bi(x, h) - ai(x, h) * pp[i]);
            qq[i + 1] = (ai(x, h) * qq[i] - di(x, h)) / (bi(x, h) - ai(x, h) * pp[i]);
        }
        y1[n] = (an * qq[n] - dn) / (bn - an * pp[n]);
        for (; i > 0; i--)
        {
            y1[i - 1] = pp[i] * y1[i] + qq[i];
        }
        n *= 2;
        h = (b - a) / n;
        b0 = -(h * a1 - a2);
        c0 = a2;
        d0 = h * A;
        an = -b2;
        bn = -(h * b1 + b2);
        dn = h * B;
        y2 = new double[n + 1];
        pp = new double[n + 1];
        qq = new double[n + 1];
        i = 0;
        pp[i] = c0 / b0;
        qq[i] = (-d0) / b0;
        for (x = a; i < n; i++, x += h)
        {
            pp[i + 1] = ci(x, h) / (bi(x, h) - ai(x, h) * pp[i]);
            qq[i + 1] = (ai(x, h) * qq[i] - di(x, h)) / (bi(x, h) - ai(x, h) * pp[i]);
        }
        y2[n] = (an * qq[n] - dn) / (bn - an * pp[n]);
        for (; i > 0; i--)
        {
            y2[i - 1] = pp[i] * y2[i] + qq[i];
        }
    } while (!checkRunge(y1, y2, n / 2 + 1));
    h = (b - a) / saverN;
    x = a;
    int j = 0;
    for (int i = 0; i <= n; i += n / saverN){
        cout << "x[" << j << "]: " << x << " \ty[" << j << "]: " << y2[i] << endl;
        j++;
        x += h;
    }

    cout << endl;
    cout << "n = " << saverN << endl;
    cout << "h = " << h << endl;

    return 0;
}