Написать программу решения системы линейных алгебраических уравнений

#include <stdio.h>
 
void create_vector(double** a, int n) { *a = new double[n]; }
void create_matrix(double*** a, int n, int m) { int i; *a = new double*[n]; for(i=0;i<n;i++) (*a)[i] = new double[m]; }
void print_matrix(double** a, int n, int m) { int i,j; for(i = 0; i < n; i++) { for(j = 0; j < m; j++) printf("%f ", a[i][j]); printf("\n"); }  }
void print_vector(double* a, int n) { int i; for(i = 0; i < n; i++) printf("%f ", a[i]); printf("\n"); }
 
void mul(double** a1, double** a2, double** r, int n, int m)
{
 int i, j, k; 
 double sum;
 for(i = 0; i < n; i++) // строка первой матрицы 
  for(j = 0; j < m; j++) // столбец второй матрицы
  { 
   sum = 0.0; 
   for(k = 0; k < m; k++) // столбец первой матрицы   
    sum += a1[i][k] * a2[k][j];      
   r[i][j] = sum;      
  }  
}
 
 
int main()
{
int i,j,k,n(3),m(3);
double sum;
double **a;
double **l;
double **u;
double **r;
double *b;
double *y;
double *x;
 
create_matrix(&a,3,3);
create_matrix(&l,3,3);
create_matrix(&u,3,3);
create_matrix(&r,3,3);
create_vector(&b,3);
create_vector(&y,3);
create_vector(&x,3);
 
a[0][0] = 2; a[0][1] = -0.24; a[0][2] = 1;
a[1][0] = 3; a[1][1] = 5; a[1][2] = -2;
a[2][0] = 1; a[2][1] = -4; a[2][2] = 10;
 
b[0] = -3; b[1] = 1; b[2] = 0;
 
printf("a\n"); print_matrix(a, 3, 3); printf("\n");
printf("b\n"); print_vector(b,3); printf("\n");
 
// habr.com/ru/sandbox/35982/
for(i=0;i<3;i++)
 for(j=0;j<3;j++)
  u[i][j] = a[i][j];
 
 
for(i=0;i<3;i++)
 for(j=i;j<3;j++)
  l[j][i] = a[i][i];
 
for(k = 1; k < 3; k++)
{
 for(i = k-1; i < 3; i++)
  for(j = i; j < 3; j++)
   l[j][i] = u[j][i]/u[i][i];
 
 for(i = k;i < 3; i++)
  for(j = k-1;j < 3; j++)
   u[i][j] = u[i][j] - l[i][k-1]*u[k-1][j];
   
}
 
printf("u\n"); print_matrix(u, 3, 3); printf("\n");
printf("l\n"); print_matrix(l, 3, 3); printf("\n");
mul(l,u,r,3,3);
printf("r\n"); print_matrix(r, 3, 3); printf("\n");
 
y[0] = b[0] / l[0][0];
y[1] = (b[1] - l[1][0]*y[0]) / l[1][1];
y[2] = (b[2] - l[2][0]*y[0] - l[2][1]*y[1]) / l[2][2];
 
printf("y\n"); print_vector(y,3); printf("\n");
 
x[2] = y[2] / u[2][2]  ;
x[1] = (y[1] - u[1][2]*x[2]) / u[1][1];
x[0] = (y[0] - u[0][2]*x[2] - u[0][1]*x[1]) / u[0][0];
 
printf("x\n"); print_vector(x,3); printf("\n");
 
printf("%f\n", a[0][0]*x[0]+a[0][1]*x[1]+a[0][2]*x[2]);
printf("%f\n", a[1][0]*x[0]+a[1][1]*x[1]+a[1][2]*x[2]);
printf("%f\n", a[2][0]*x[0]+a[2][1]*x[1]+a[2][2]*x[2]);
 
return 0;
}

y[0] = b[0] / l[0][0];
    y[1] = (b[1] - l[1][0] * y[0]) / l[1][1];
    y[2] = (b[2] - l[2][0] * y[0] - l[2][1] * y[1]) / l[2][2];
    y[3] = (b[3] - l[3][0] * y[0] - l[3][1] * y[1]) / l[3][3];
 
    printf("y\n"); print_vector(y, 4); printf("\n");
 
//Вот тут завис...
    x[2] = y[2] / u[2][2]; 
    x[1] = (y[1] - u[1][2] * x[2]) / u[1][1];
    x[0] = (y[0] - u[0][2] * x[2] - u[0][1] * x[1]) / u[0][0];
 
    printf("x\n"); print_vector(x, 4); printf("\n");
 
    printf("%f\n", a[0][0] * x[0] + a[0][1] * x[1] + a[0][2] * x[2] + a[0][3] * x[3]);
    printf("%f\n", a[1][0] * x[0] + a[1][1] * x[1] + a[1][2] * x[2] + a[1][3] * x[3]);
    printf("%f\n", a[2][0] * x[0] + a[2][1] * x[1] + a[2][2] * x[2] + a[2][3] * x[3]);
    printf("%f\n", a[3][0] * x[0] + a[3][1] * x[1] + a[3][2] * x[2] + a[3][3] * x[3]);
 
    std::cout << "Press any key." << "\n";
    std::cin.get();
    return 0;