Запутался с методом Рунге-Кутты

[Morpho]

Всем привет
Имеются вот такие уравнения:


Мне требуется проинтегрировать их методом Рунге-Кутты 4 порядка.
Написал вот такой код.

for q := 1 to vols do
    begin
      h0:=h+step;
      ex0:=ex+step;
      ey0:=ey+step;
      G0:=G+step;

      k11 := step * (h0*h0/ksi * F * sin(lambda));
      k21 := step * ((h0+0.5*step)/ksi * F * sin(lambda));
      k31 := step * ((h0+0.5*step)/ksi * F * sin(lambda));
      k41 := step * ((h0+step)/ksi * F * sin(lambda));

      k12 := step * (h*F*( sin(G0)*cos(lambda) + ( cos(G0)*(ksi-1) + ex0 ) * sin(lambda) ));
      k22 := step * (h*F*( sin(G0)*cos(lambda) + ( cos(G0)*(ksi-1) + (ex0+0.5*k12) ) * sin(lambda) ));
      k32 := step * (h*F*( sin(G0)*cos(lambda) + ( cos(G0)*(ksi-1) + (ex0+0.5*k22) ) * sin(lambda) ));
      k42 := step * (h*F*( sin(G0)*cos(lambda) + ( cos(G0)*(ksi-1) + (ex0+k32) ) * sin(lambda) ));

      k13 := step * (h*F*(-cos(G0)*cos(lambda) + ( sin(G0)*(ksi+1) + ey0 ) * sin(lambda)) );
      k23 := step * (h*F*(-cos(G0)*cos(lambda) + ( sin(G0)*(ksi+1) + (ey0+0.5*k13) * sin(lambda)) ));
      k33 := step * (h*F*(-cos(G0)*cos(lambda) + ( sin(G0)*(ksi+1) + (ey0+0.5*k23) * sin(lambda)) ));
      k43 := step * (h*F*(-cos(G0)*cos(lambda) + ( sin(G0)*(ksi+1) + (ey0+k33) * sin(lambda)) ));

      k14 := step * ( sqr(1 + ex0*cos(G0) + ey0*sin(G0)) / (h*h*h * mu) );
      k24 := step * ( sqr(1 + ex0*cos(G0+0.5*k14) + ey0*sin(G0+0.5*k14)) / (h*h*h * mu) );
      k34 := step * ( sqr(1 + ex0*cos(G0+0.5*k24) + ey0*sin(G0+0.5*k24)) / (h*h*h * mu) );
      k44 := step * ( sqr(1 + ex0*cos(G0+k34) + ey0*sin(G0+k34)) / (h*h*h * mu) );

      h := h0 + (k11 + 2*k21 + 2*k31 + k41) / 6;
      ex := ex0 + (k12 + 2*k22 + 2*k32 + k42) / 6;
      ey := ey0 + (k13 + 2*k23 + 2*k33 + k43) / 6;
      G := G0 + (k14 + 2*k24 + 2*k34 + k44) / 6;
end;

И вот вопрос: в методе Рунге-Кутты при расчёте коэффициентов К2, К3 и К4 требуется увеличивать аргумент на пол шага и на шаг... В этом месте я ошибся, видимо... Можете подсказать, как исправить ситуацию?

[Aspirisha]

Да, вроде все достаточно прозрачно. Нужны разные интервалы интегрирования. Разумно принять h например за независимую переменную, ее интегрировать с шагом step / 4. G - с шагом step/2, ex и ey с шагом step. Соответственно размеры массивов будут разные, и в целом цикл for надо подтесать (например, в лоб if step mod 4 == 0...)

[Morpho]

А вот ещё такой вопрос. Во всех методах (Эйлера и Рунге-Кутты 4 порядка) в общем виде дифур выглядит как dy/dt = f(x, y). Если с пониманием, что при программировании в моём случае брать вместо у вопросов нет (брать ту переменную, которая стоит слева), то что принимать за икс? Все остальные переменные? И их все увеличивать на шаг? (в методе Эйлера)

[Aspirisha]

Ну в общем да. У вас система дифуров:
dh/dt = F1(h, ex, ey, G, t)
dex/dt = F2(h, ex, ey, G, t)
...
X - видимо считается при такой записи как у вас за вектор всех переменных кроме той что слева в числителе.