#include <iostream>
#include <cmath>

using namespace std;

struct Point {
    double x, y;

    Point() {x = y = 0;}
    Point(double _x, double _y) : x(_x), y(_y) {};
};

double checkEps(Point* s1, Point* s2, int n) {
    double eps = -1, curEps;

    for (int i = 0; i < n; i++) {
        curEps = fabs(s1[i].y - s2[i * 2].y) / 15;
        eps = max(eps, curEps);
    }
    return eps;
}

double pi(double x) {return x;}

double qi(double x) {return 1;}

double fi(double x) {return x + 1;}

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

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

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

double di(double x, double h) {return h * h * fi(x);}

Point* createPoint(double xa, double xb, double a1, double a2, double a, double b1, double b2, double b, int n) {
    Point* res = new Point[n + 1];
    double pp[n + 1], qq[n + 1], h = (xb - xa) / n, b0 = -(h * a1 - a2), c0 = a2, d0 = h * a, an = -b2, bn = -(h * b1 + b2), dn = h * b;
    int i = 0;

    pp[0] = c0 / b0; qq[0] = -d0 / b0;

    for (double x = xa; 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]);
    }
    res[n] = Point(xb, (an * qq[n] - dn) / (bn - an * pp[n]));

    for (i; i > 0; i--)
        res[i - 1] = Point(xb - h * (n - i + 1), pp[i] * res[i].y + qq[i]);

    return res;
}

void solve(double xa, double xb, double a1, double a2, double a, double b1, double b2, double b, double eps) {
    Point *s1, *s2;
    int n = 10;

    do {
        s1 = createPoint(xa, xb, a1, a2, a, b1, b2, b, n);
        n *= 2;
        s2 = createPoint(xa, xb, a1, a2, a, b1, b2, b, n);
    } while (checkEps(s1, s2, n / 2) > eps);

    cout << "n = " << n << endl;
    cout << endl;
    cout.precision(5);
    for (int i = 0; i <= n; i++) {
        cout << "x[" << i << "] = " << s2[i].x << "\ty[" << i << "] = " << s2[i].y << endl;
    }
}

int main() {
    solve(0.5, 0.8, 1., 2., 1., 0., 1., 1.2, 1e-4);

    return 0;
}