Квадратичная аппроксимация, ошибка при решении СЛАУ методом Гаусса

// функция решения слау методом Гаусса, а - матрица, х - куда пихать результат, н - порядок матрицы
// у нас везде нумерация начинается с 1 а не 0, тоисть для матрици 3х4(для нашой квадратичной 
// аппроксимации) матрица а будет 4х5, массив х будет из 4-х ячеек, н = 4), нулевые ячейки, рядки и .
//столбци не учитываются
 
void Gays( double **a, double *x, int n )
{
  int m = n+1;
  int i,j;
  for( i = 1 ; i < n-1 ; i++ )
    for( j = i ; j < m-1 ; j++ ){
      int t = i;
      while( t < n-1 )
      {
        a[t+1][j+1] -=a [t+1][i] / a[i][i] * a[i][j+1];
        t++;
      }
   }
 
   for(int i = n-1 ; i >= 1 ; i-- )
   {
     x[i] = a[i][m-1] / a[i][i];
     for( int j = m-2 ; j >= i && i > 1 ; j-- )
       a[i-1][m-1] -= a[i-1][j] * x[j];
   }
}
//-----------------------------------------------------------------------
// функция для аппроксимации
double fx2( double x ){ return( 0.0 + 0.0*x + pow(x, 0.25) ); }
//  функция для получения значения У при известных коефициентах(полученых полсле решения СЛАУ Гауссом) матрици
double fx2( double x, double a1, double a2, double a3 ){ return( a1 + a2*x + a3*x*x ); }
//---------------------------------------------------------------------------
 
void f()
{
  int n = 4;                               // порядок матрици+1 
  int m = n+1;
 
  double *k  = new double[m];
  double **a = new double*[n];
  for(int i=0;i<n;i++)
    a[i] = new double[m];
 
  double Xp = StrToFloat(  Edit2->Text ); //  Хр = 1
  double Xk = StrToFloat(  Edit3->Text ); //  Хк = 10
  double shag= StrToFloat( Edit4->Text );//  = 0.5
  int iter = (Xk-Xp)/shag;
 
  double *x = new double[ iter ];
  double *y = new double[ iter ];
 
 
  int moo = 0;
  for( double i = Xp ; i < Xk ; i+=shag )  //  получаем експерементальные данные( наша функция + шум )
  {
    x[moo] = i;
    y[moo] = fx2(i) + random(200)*0.001 - 0.1 ;
    moo++;
  }
    //-------------------
  double xi  = 0.0;
  double x2i = 0.0;
  double yi  = 0.0;
  double x3i = 0.0;
  double xyi = 0.0;
  double x4i = 0.0;
  double xy2i= 0.0;
 
 
  for( int i = 0 ; i < iter ; i++ ) // вычисляем данные для формирования матрици 
  {
    xi  += x[i];
    x2i += x[i]*x[i];
    yi  += y[i];
    x3i += x[i]*x[i]*x[i];
    xyi += x[i]*y[i];
    x4i += x[i]*x[i]*x[i]*x[i];
    xy2i += x[i]*y[i]*y[i];
  }
    
  //  заполняем матрицу как на картинке
  a[1][1] = iter;
  a[1][2] = xi;
  a[1][3] = x2i;
  a[1][4] = yi;
  a[2][1] = xi;
  a[2][2] = x2i;
  a[2][3] = x3i;
  a[2][4] = xyi;
  a[3][1] = x2i;
  a[3][2] = x3i;
  a[3][3] = x4i;
  a[3][4] = xy2i;
   //-------------------
 
 Gays( a, k, n );  // решаем СЛАУ(созданную матрицу) и резульрат пихаем в к
 
 for(int i = 1 ; i < n ; i++ )
   Memo1->Lines->Add( FloatToStr(k[i]) );  //печатаем полученные результаты
 
 Chart1->Series[1]->Clear();
 Chart1->Series[0]->Clear();
 Chart1->Series[2]->Clear();
 for( int i = 0 ; i < iter ; i++ )
    Chart1->Series[0]->AddXY(x[i],y[i]);  // рисуем экспериментальные данные
 
 for( double i = Xp ; i < Xk ; i += shag )
 {
    Chart1->Series[1]->AddXY(i, fx2(i) );                         // рисуем настоящее значение функции
    Chart1->Series[2]->AddXY(i, fx2(i,k[1], k[2], k[3] ));    // рисуем аппроксимуючу   функцию 
 }
 // в результате либерда, скрин с результатом ниже
 
     //-------------------
  delete []k;
  for(int i=0;i<n;i++)
    delete []a[i];
  delete []a;
 
}
//---------------------------------------------------------------------------