Метод Милна

[PhotOn]

Приветствую вас. Очень помощь нужна мне © Йода XD. Нужно модернизировать программку, чтобы можно было производить расчет по нескольким уравнениям, а не по одному, заданному программой. Ниже вложено изображение с примерами. Все что нужно сделать в программе: сделать возможным ввод коэффициентов уравнений и выбора степени х и у.

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));
}
}

//------