Метод Милна. Расчет по нескольким уравнениям.

float TForm1::f1(float xa, float ya, float yb)
{
return 2*xa*xa+2*ya+yb;
}
// Второе уравнение системы
float TForm1::f2(float xa, float ya, float yb)
{
return 1-2*xa*xa+2*ya-yb;
}
void __fastcall TForm1::Button1Click(TObject *Sender)
{
float h;   // Шаг
float a,b,k1,k2,k3,k4;
float r1,r2,r3,r4;
float eps,abs_pogr;   // eps - точность
float zkor[11],zpr[11],ypr[11],ykor[11],x[11],y1[11],y2[11];
int i=0;
// Ввод концов отрезка
a=StrToFloat(Edit1->Text);
b=StrToFloat(Edit2->Text);
x[0]=a;
// Начальные условия:
y1[0]=StrToFloat(Edit3->Text);
y2[0]=StrToFloat(Edit4->Text);
// Шаг
h=StrToFloat(Edit5->Text);
// Точность
eps=StrToFloat(Edit6->Text);
// Решение системы уравнений методом Рунге-Кутта
while (i<=3)
{
k1=h*f1(x[i],y1[i],y2[i]);
r1=h*f2(x[i],y1[i],y2[i]);
k2=h*f1(x[i]+h/2,y1[i]+k1/2,y2[i]+r1/2);
r2=h*f2(x[i]+h/2,y1[i]+k1/2,y2[i]+r1/2);
k3=h*f1(x[i]+h/2,y1[i]+k2/2,y2[i]+r2/2);
r3=h*f2(x[i]+h/2,y1[i]+k2/2,y2[i]+r2/2);
k4=h*f1(x[i]+h,y1[i]+k3,y2[i]+r3);
r4=h*f2(x[i]+h,y1[i]+k3,y2[i]+r3);
 
y1[i+1]=y1[i]+(k1+2*k2+2*k3+k4)/6;
y2[i+1]=y2[i]+(r1+2*r2+2*r3+r4)/6;
 
x[i+1]=x[i]+h;
i=i+1;
}
i=4;
// Решение системы уравнений методом Милна
while (x[i]<=b+h)
{
// Для предсказания используется первая формула Милна
ypr[i]=y1[i-4]+(4*h)/3*(2*f1(x[i-3],y1[i-3],y2[i-3])-f1(x[i-2],y1[i-2],y2[i-2])+2*f1(x[i-1],y1[i-1],y2[i-1]));
// Для уточнения - вторая формула Милна
ykor[i]=y1[i-2]+(h/3)*(f1(x[i-2],y1[i-2],y2[i-2])+4*f1(x[i-1],y1[i-1],y2[i-1])+f1(x[i],ypr[i],y2[i]));
// Для второго уравнения
zpr[i]=y2[i-4]+(4*h)/3*(2*f2(x[i-3],y1[i-3],y2[i-3])-f2(x[i-2],y1[i-2],y2[i-2])+2*f2(x[i-1],y1[i-1],y2[i-1]));
zkor[i]=y2[i-2]+(h/3)*(f2(x[i-2],y1[i-2],y2[i-2])+4*f2(x[i-1],y1[i-1],y2[i-1])+f2(x[i],zpr[i],y2[i]));
 
abs_pogr=abs(ykor[i]-ypr[i])/29;
if (abs_pogr>eps) y1[i]=ykor[i];
else y1[i]=ypr[i];
// Абсолютная погрешность
abs_pogr=abs(zkor[i]-zpr[i])/29;
// Контроль точности полученного результата
if (abs_pogr>eps) y2[i]=zkor[i];
else y2[i]=zpr[i];
 
x[i+1]=x[i]+h;
i=i+1;
}
// Вывод решения системы уравнений в виде функций y1[i], y2[i]
ListBox1->Items->Add("x[i]         y1[i]            y2[i]");
for (i=0;i<=10;i++)
{
ListBox1->Items->Add(FloatToStrF(x[i],ffFixed,3,1)+"      "+
FloatToStrF(y1[i],ffFixed,10,6)+
"      "+FloatToStrF(y2[i],ffFixed,10,6));
}
}
 
//------