Метод рунге-кутты-мерсона - C (СИ)

Формулировка задачи:

нужно написать программу для решения систем из n дифф уравнений методом рунге-кутты-мерсона(он же пятиэтапный метод рунге кутта 4 порядка). вот кусок кода,который считает для одного уравнения.написан не идеально,только учусь..вот не могу понять как сделать для системы..уже сижу дня три и уже голову сломал..есть еще куски,пробовал написать для двух уравнений,но такое не принимают... спасите меня...хоть какой-нибудь идейкой...просто исчерпал все свои знания. у нас есть система вида y1'=a1*y1+b1*y2+..+N1*yn y2'=a2*y1+b2*y2+..+N2*yn ..................................... yn'=an*y1+bn*y2+..+Nn*yn и есть начальные условия,как для задачи коши

do
    {
        k1=f(x[i],y[i]);
        k2=f(x[i]+h/3,y[i]+h/3*k1);
        k3=f(x[i]+h/3,y[i]+h/6*k1+h/6*k2);
        k4=f(x[i]+h/2,y[i]+h/8*k1+3*h/8*k2);
        k5=f(x[i]+h,y[i]+h/2*k1-3*h/2*k3+2*h*k4);
        yx=y[i]+h/2*(k1-3*k3+4*k4);
        y=y[i]+h/6*(k1+4*k4+k5);
        r=0.2*abs(y-yx);
        if(r>e) h=h/2;
        else 
        {
            if((r<=e/64)&&(2*h<h0)) h=2*h;
            else h=h;
        }
        while(r>e)
        {
            k1=f(x[i],y[i]);
            k2=f(x[i]+h/3,y[i]+h/3*k1);
            k3=f(x[i]+h/3,y[i]+h/6*k1+h/6*k2);
            k4=f(xi]+h/2,y[i]+h/8*k1+3*h/8*k2);
            k5=f(x[i]+h,y[i]+h/2*k1-3*h/2*k3+2*h*k4);
            yx=y[i]+h/2*(k1-3*k3+4*k4);
            y=y[i]+h/6*(k1+4*k4+k5);
            r=0.2*abs(y-yx);
            if(r>e) h=h/2;
            else 
            {
                if((r<=e/64)&&(2*h<h0)) h=2*h;
                else h=h;
            }
        }
 
        y[i+1]=y;
        i=i+1;
        x[i]=x[i-1]+h;
        
    }
    while (x[i-1]<=xn);

Решение задачи: «Метод рунге-кутты-мерсона»

Листинг программы

double X,X9,H, EH,EPS,H0,P,D,X1,Y[8];
//------------------------------------------------------ ГЇГ°Г®ГЁГ§ГўГ®Г¤Г*Г»ГҐ *>
void RP(double X,double *Y,double *F)
{
   F[0] = Y[1];
   F[1] = P * (1.0 - Y[0]*Y[0])*Y[1]-Y[0];
}
//------------------------------------------------ ГЊГ…Г’ГЋГ„ ГђГ“ГЌГѓГ…-ГЉГ“Г’Г’Г›-ГЊГ…ГђГ‘ГЋГЌГЂ *>
void RKM(int N,double *X,double *H,double *E, double *Y)
{
        int I;
        double H3,H4,R = 0.0,A,Z[8],F0[8],F[8],K1[8],K3[8];
        RP(*X,&Y[0],&F[0]);
        for (I=0;I<N ;I++)
           Z[I] = Y[I];
        do
        {
                H3 = *H / 3.0;
                H4 = 4 * H3;
                for (I=0; I<N; I++)
                {
                        K1[I] = H3 * F0[I];
                        Y[I] = Z[I] + K1[I];
                }
                RP(*X + H3,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                   Y[I] = Z[I] + (K1[I] + H3*F[I])/2;
                RP(*X + H3,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K3[I] = *H * F[I];
                        Y[I] = Z[I] + 0.375 * (K1[I] + K3[I]);
                }
                RP((*X + *H/2),&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K1[I] += H4 * F[I];
                        Y[I] = Z[I] + 1.5 * (K1[I] + K3[I]);
                }
                RP(*X + *H/2,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K1[I] += H4 * F[I];
                        Y[I] = Z[I] + 1.5 * (K1[I] - K3[I]);
                }
                RP(*X + *H,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        A = H3 * F[I];
                        Y[I] = Z[I] + (K1[I] + A)/2;
                        A = 2 * K1[I] - 3.0 * K3[I] - A;
                        if ( Y[I] != 0.0)
                                A /= Y[I];
                        if ( abs(A) > R )
                                R=abs(A);
                }
                *H = *H/2;
        }
        while(R<*E);
        *H *= 2;
        *X += *H;
        if (32*R < *E)
                *H *= 2;
}
void __fastcall TForm1::run() 
{
//входящие данные
        X=Edit1->Text.ToDouble();
        X9=Edit2->Text.ToDouble();
        H=Edit3->Text.ToDouble();
        EH=Edit4->Text.ToDouble();
        EPS=Edit5->Text.ToDouble();
        H0=Edit6->Text.ToDouble();
        Y[0]=Edit7->Text.ToDouble();
        Y[1]=Edit8->Text.ToDouble();
        P=Edit9->Text.ToDouble();
//входящие данные
        X1 = X;
        RKM(2,&X,&H0,&EPS,Y);
        H0 =  Sign (H);
        D  = X1 + H - X;
        if (D < 0.0 == (H>0.0))
                H0 = D;
}

Объяснение кода листинга программы

Входные данные:
- X = Edit1->Text.ToDouble()
- X9 = Edit2->Text.ToDouble()
- H = Edit3->Text.ToDouble()
- EH = Edit4->Text.ToDouble()
- EPS = Edit5->Text.ToDouble()
- H0 = Edit6->Text.ToDouble()
- Y[0] = Edit7->Text.ToDouble()
- Y[1] = Edit8->Text.ToDouble()
- P = Edit9->Text.ToDouble()
Код реализует метод RKM (Runge-Kutta-Merson) для решения обыкновенного дифференциального уравнения.
Переменные, используемые в методе RKM:
- Z[i] - временная переменная для хранения промежуточных значений
- F0[i] - массив для хранения начальных значений функций
- F[i] - массив для хранения текущих значений функций
- K1[i] - массив для хранения первых приближений
- K3[i] - массив для хранения третьих приближений
- A - временная переменная для хранения значения величины приращения
Выходные данные:
- H0 - полученное приближение решения
- D - ошибка при решении
В методе TForm1::run() происходит следующее:
- Входные данные считываются из текстовых полей и преобразуются в числовые значения.
- Вызывается функция RKM с параметрами N=2 и указанием входных переменных X, H0, EPS, Y.
- Значение H0 изменяется на знак числа H.
- Вычисляется значение D как разность между X1 и X, где X1 - предыдущее значение X.
- Если D меньше нуля и H больше нуля, то H0 присваивается значение D.