include iostream include cmath include vector using namespace std type

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


using namespace std;

typedef std::vector<double> vec;
typedef std::vector<vec> Matrix;
const double eps = 1e-5;

const double a = 2.7;
const double b = 3.2;
const double alpha0 = 2;
const double alpha1 = -0.2;
const double A = 3;
const double beta0 = -7;
const double beta1 = 0.4;
const double B = 6;

const double h = 0.05;
int n = int((b - a) / h);

Matrix S(n + 1, vec(n + 2, 0));
vec y_i(n + 1, 0);
double c1;
double c2;

double p(double x)
{
    return -6.0;
}

double q(double x)
{
    return 7.0;
}

double f(double x)
{
    return 2*x-1;
}

double x_i(int i)
{
    return a + i * h;
}

double m_i(int i)
{
    return -2.0 + h * p(x_i(i));
}

double n_i(int i)
{
    return 1.0 - h * p(x_i(i)) + h * h * q(x_i(i));
}

void solve()
{
    for (int i = 0; i < n + 1; i++)
    for (int j = 0; j < n + 2; j++)
        S[i][j] = 0.0;

    for (int i = 0; i < n - 1; i++)
    {
        S[i][i + 2] = 1.0;
        S[i][i + 1] = m_i(i);
        S[i][i] = n_i(i);
        S[i][n + 1] = f(x_i(i)) * h * h;
    }

    S[n - 1][0] = alpha0 - alpha1 / h;
    S[n - 1][1] = alpha1 / h;
    S[n - 1][n + 1] = A;

    S[n][n] = beta0 + beta1 / h;
    S[n][n - 1] = -beta1 / h;
    S[n][n + 1] = B;
}


void triadiag()
{
    for (int i = 0; i < n + 1; i++)
    {
        int u = i;
        while (u < n + 1 && fabs(S[u][i]) < eps)
            u++;

        if (u != i && u < n + 1 && fabs(S[u][i]) >= eps)
        for (int j = 0; j < n + 2; j++)
        {
            double tmp = S[i][j];
            S[i][j] = S[u][j];
            S[u][j] = tmp;
        }

        if (fabs(S[i][i]) >= eps)
        {
            double rem = S[i][i];
            for (int j = i; j < n + 2; j++)
                S[i][j] /= rem;

            for (int v = i + 1; v < n + 1; v++)
            {
                double mul = S[v][i];
                for (int j = i; j < n + 2; j++)
                    S[v][j] -= S[i][j] * mul;
            }
        }

    }

    for (int i = n; i >= 0; i--)
    {
        y_i[i] = S[i][n + 1];
        for (int j = n; j > i; j--)
            y_i[i] -= S[i][j] * y_i[j];
    }
}

int main()
{
    solve();
    triadiag();
    for (int i = 0; i < n + 1; i++)
    {
        cout << "y" <<'(' << x_i(i) << ") = " << y_i[i] << endl;
    }

}