diff equation accepted

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

using namespace std;

#define xtype double

#define N 11

//y''(x) + p(x)y'(x) + q(x)y(x) = f(x)
//y(a) = A
//y(b) = B

xtype y(xtype x) { return 3*x*x*x; }
xtype dy(xtype x) { return 9*x*x; }
xtype ddy(xtype x) { return 18*x; }

xtype pFunc(xtype x) { return 1; }
xtype qFunc(xtype x) { return -1; }


xtype f[N], p[N], q[N];

xtype x_A = 1;
xtype x_B = 2;

xtype y_A = y(x_A);
xtype y_B = y(x_B);

xtype h = (x_B - x_A) / (N-1);



void solveMatrix(int n, xtype* a, xtype* b, xtype* c, xtype* d, xtype* x) {

    int i;
    xtype alpha[n], beta[n];
    alpha[0] = -c[0] / b[0];
    beta[0] = d[0] / b[0];

    for (i = 1; i < n; i++) {
        alpha[i] = (-c[i]) / (a[i] * alpha[i-1] + b[i]);
        beta[i] = (d[i] - a[i] * beta[i-1]) / (a[i] * alpha[i-1] + b[i]);
    }

    x[n-1] = beta[n-1];

    for (i = n-2; i >= 0; i--) {
        x[i] = alpha[i] * x[i+1] + beta[i];
    }
}

void countFunc(int n) {
    for (int i = 0; i < n; i++) {

        xtype x = x_A + i*h;

        f[i] = ddy(x) + pFunc(x)*dy(x) + qFunc(x)*y(x);
        p[i] = pFunc(x);
        q[i] = qFunc(x);
    }

}

void fillAndSolveMatrix(int n) {

    xtype a[n], b[n], c[n], d[n], e[n], res[n];

    /*    a[0] = 0;
    b[0] = q[1]*h*h - 2;
    c[0] = 1 + p[1]*h/2;
    d[0] = f[1]*h*h - y_A*(1 - p[1]*h/2);

    a[n-3] = 1 - p[n-2]*h/2;
    b[n-3] = q[n-2]*h*h - 2;
    c[n-3] = 0;
    d[n-3] = f[n-2]*h*h - y_B*(1 + p[n-2]*h/2);*/


    for (int i = 0; i < n; i++) {

        a[i] = 1 - p[i+1]*h/2;
        b[i] = q[i+1]*h*h - 2;
        c[i] = 1 + p[i+1]*h/2;
        d[i] = f[i+1]*h*h;
    }

    a[0] = 0;
    b[0] = q[1]*h*h - 2;
    c[0] = 1 + p[1]*h/2;
    d[0] = f[1]*h*h - (1 - p[1]*h/2)*y_A;

    a[n-2] = 1 - p[n-1]*h/2;
    b[n-2] = q[n-2]*h*h - 2;
    c[n-2] = 0;
    d[n-2] = f[n-1]*h*h - (1 + p[n-1]*h/2)*y_B;



    solveMatrix(n-1, a, b, c, d, e);

    for (int i = 0; i < n-1; i++) res[i+1] = e[i];
    res[0] = y_A;
    res[n] = y_B;

    for (int i = 0; i <= n; i++) {
        printf("x[%d]=%.6f y[%d]=%.6f res[%d]=%.6f d[%d]=%.6f\n",
                i, x_A+i*h, i, y(x_A+i*h), i, res[i], i, fabs(y(x_A+i*h)-res[i]));
    }
}



int main()
{

    countFunc(N);
    fillAndSolveMatrix(N-1);

    return 0;
}