Ошибка в методе Рунге-Кутты

[illuminates]

Здравствуйте! Помогите пожалуйста найти неполадки в программе. Дело в том, что программа должна решать диффуру методом Рунге-Кутты, точнее, заполнять три динамических массива х - переменная, у - искомая функция, z - первая производная функции у. На выходе она заполняет файл данными этих массивов. Но значение у одинаково во всех строках. При проверке оказалось, что коэффициенты k0, k1, k2, k3 на каждом шаге равны нулю! Хотя коэффициенты l0, l1, l2, l3 и шаг h нулю не равны.
Программа вроде-бы прозрачная, но ошибку найти не получается. К слову, подпрограмма метода Рунге-Кутты хорошо работала в другой моей задаче давала верные решения.

int main (void)
{
   double alf_0, alf_1 = 0.9, alf_2 = 1.05, epsilon = 0.001, x_0 = 0, x_n = 1, z_0 = 0, h, A, B, C;
   short n = 20;
   double *x = new double [n+1];
   double *y = new double [n+1];
   double *z = new double [n+1];
   short i;
   
   h = (x_n - x_0)/ n;
   cout << h<<"\n";
   /*do
   {*/
      alf_0 = (alf_1 + alf_2)/ 2;
      
      Run_Kut (x, y, z, x_0, alf_1, z_0, h, n);
      A = y[n] + 2 * z[n];
         
      for (i = 0; i<n; i++)
         cout << x[i] <<"\t"<< y[i] <<"\t"<< z[i] <<"\n";

       Run_Kut (x, y, z, x_0, alf_0, z_0, h, n); 
      C = y[n] + 2 * z[n];

      if ((A * C)> 0) alf_1 = alf_0;
      else alf_2 = alf_0;
   /*}
   while (alf_2 - alf_1 > epsilon);
   
   ofstream out;
    out.open ("D:\\Задача 6_1.txt");*/
   for (i = 0; i<n; i++)
      cout << x[i] <<"\t"<< y[i] <<"\t"<< z[i] <<"\n";
      /*out << x[i] <<"\t"<< y[i] <<"\t"<< z[i] <<"\n";
   out.close();*/
   delete[] x, y, z;
   getch();
   return 0;
}
/**/
void Run_Kut (double *X, double *Y, double *Z, double x_0, double y_0, double z_0, double h, short n)
{
   double k0, k1, k2, k3,   l0, l1, l2, l3;
   int i;
   //начальные условия
   X[0] = x_0; Y[0] = y_0; Z[0] = z_0;
   for (i = 1; i<= n; i++)
   {
      k0 = Fy (Z[i-1], 0, 0);                   l0 = Fz (X[i-1],       Y[i-1],             Z[i-1]);
      k1 = Fy (Z[i-1] + h * l0/ 2, 0, 0);  l1 = Fz (X[i-1] + h/2, Y[i-1] + h * k0/ 2, Z[i-1] + h * l0/ 2);
      k2 = Fy (Z[i-1] + h * l1/ 2, 0, 0);  l2 = Fz (X[i-1] + h/2, Y[i-1] + h * k1/ 2, Z[i-1] + h * l1/ 2);
      k3 = Fy (Z[i-1] + h * l2, 0, 0);      l3 = Fz (X[i-1] + h,   Y[i-1] + h * k2,    Z[i-1] + h * l2);
      // для проверки коэффициентов
      /**/
      cout << k0 <<"   "<< k1 <<"   "<< k2 <<"   "<< k3 <<"\n";
      cout << l0 <<"   "<< l1 <<"   "<< l2 <<"   "<< l3 <<"\n";

      X[i] = X[i-1] + h;
      Y[i] = Y[i-1] + (h/ 6) * (k0 + 2 * k1 + 2 * k2 + k3);
      Z[i] = Z[i-1] + (h/ 6) * (l0 + 2 * l1 + 2 * l2 + l3);
   }
}

double Fz (double xi, double yi, double zi)
{
   return ( -2 * xi * zi + xi * yi + exp(-xi) )/(pow(xi, 2) + 1);
}
double Fy (double xi, double yi, double zi)
{
   return zi;
}

Заранее спасибо!

[ewert]

Вы просто перепутали порядок аргументов при обращении к функции Fy.

[illuminates]

Спасибо!