Решение краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей

#include <cstdlib> 
#include <iostream> 
#include <string>
using namespace std;
 
double c1, c2, h, tmp;
int n = 100, i, j, l, k;
double a[100][100];
double b[100], y[100], x[100];
double xbeg, xend;
void gauss();
double f(double);
double p(double);
double q(double);
int main()
{
    setlocale(0, "");
    n=11;
 
    xbeg=-1;
    xend=1;
    h = 0.2;
 
    for (i = 0; i < n; i++)
    {
        b[i] = 0;
        y[i] = 0;
        x[i] = xbeg + i * h;
    }
 
    //Инициализация двумерного динамического массива
    for (i = 0; i < 100; i++)
        for (j = 0; j < i + 1; j++) {
            a[i][j] = 0;
        }
    //Граничные условия
    c1 = -3;
    c2 = 2.5;
    // Правые части системы
    for (i = 0; i < n - 2; i++)
        b[i] = f(x[i]);
    //Граничные условия
    b[n - 2] = c1;
    b[n - 1] = c2;
    // Коэффициенты системы
    for (i = 0; i <= n - 3; i++)
    {
        a[i][i] = 1 / (h * h);
        a[i][i + 1] = -(2 / (h * h) + 2 * p(x[i]) / h - q(x[i]));
        a[i][i + 2] = 1 / (h * h) + 2 * p(x[i]) / h;
    }
    //Граничный условия
    a[n - 2][0] = 1;
    a[n - 1][n - 1] = 1;
    gauss();
    for (i = 0; i < n; i++)
        cout << "x(" << i << ")= " << x[i] << "  " << "y(" << i << ")=" << y[i] << endl;
    //Проверка решения на граничные условия
    cout << "Проверка решения на граничные условия" << endl;
    cout << y[0] << " " << y[n - 1] << endl;
    system("pause"); 
}
void gauss()
{
    //Прямой проход
    for (k = 0; k <= n - 2; k++)
    {
        for (i = k + 1; i <= n - 1; i++)
        {
            tmp = a[i][k];
            b[i] = b[i] - b[k] * tmp / a[k][k];
            for (j = 0; j <= n - 1; j++)
                a[i][j] = a[i][j] - a[k][j] * tmp / a[k][k];
        }
    }
    //Обратный проход
    y[n - 1] = b[n - 1] / a[n - 1][n - 1];
    for (i = n - 2; i >= 0; i = i - 1)
    {
        y[i] = b[i];
        for (j = i + 1; j <= n - 1; j++)
            y[i] = y[i] - a[i][j] * y[j];
        y[i] = y[i] / a[i][i];
    }
}
 
double f(double x)
{
    x = 2.5 * pow(x, 2) + (-3)*x;
    return x;
}
 
double p(double x)
{
    x = 2.5 * pow(x, 3)+ 3;
    return x;
}
 
double q(double x)
{
    x = (-3)*x + 2.5;
    return x;
}