float* solve_matrix(float* a, float* b, float* c, float *d) {

    float *diag_a = a;
    float *diag_b = b;
    float *diag_c = c;
    float *free_d = d;

    auto *alpha = new float[n-1];
    auto *beta = new float[n-1];
    auto *x = new float[n];


    if (diag_a[0] != 0) {
        alpha[0] = - diag_b[0]/diag_a[0];
    } else {
        std::cout << "a[0] = 0 error";
        return nullptr;
    }
    beta[0] = free_d[0] / diag_a[0];

    for (int i=1; i<n-1; i++) {
        alpha[i] = -(diag_b[i])/(alpha[i-1]*diag_c[i-1] + diag_a[i]);
        beta[i] = (free_d[i]-diag_c[i-1]*beta[i-1])/(alpha[i-1]*diag_c[i-1] + diag_a[i]);
    }

    beta[n-1] = (free_d[n-1]-diag_c[n-2]*beta[n-2])/(alpha[n-2]*diag_c[n-2] + diag_a[n-1]);

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

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

//    for (int i=0; i<n; i++) {
////        std::cout << "x[" << i << "]= " << x[i] << std::endl;
//    }
    return x;
}