Метод Рунге-Кутта 4 порядка точности

{nmax максимальное число возможных уравнений
x независимая переменная
y вектор зависимых переменных (искомое решение)
f вектор правых частей
a начальное значение зависимой переменной
b конечное значение независимой переменной
h шаг интегрирования = (b-a)/10
}
 
 
const
 
 
  nmax=8;
type
  vec=array[1..nmax] of real;
var
h,a,b,:real;
 
procedure der(x:real;y:vec; var f:vec);
begin
{вектор правых частей};
end;
procedure rk4(n: integer; x, h: real; var y: vec); {готовая процедура метода Рунге-Кутта}
var   i, j: integer;  
h1, h2, q: real;  
y0, y1, f: vec; 
begin  h1 := 0; 
h2 := h1/2;  
for i := 1 to n do begin   
y0[i] := y[i]; 
y1[i] := y[i];  
end;  
for j := 1 to 4 do begin    
der(x+h1,y,f);   
if j = 3 then h1 := h else h1 := h2;   
for i := 1 to n do begin    
q := h1*f[i]; 
y[i] := y0[i] + q;    
if j = 2 then q := 2*q;    
y1[i] := y1[i] + q/3;   
end ;  
end;  
for i := 1 to n do   y[i] := y1[i]; 
end;
 
Процесс заканчивается по достижению верхней границы например:
while x<b+0.00001 do begin
rk4(1,x,h,y);
x:=x+h;
end;

program RungeKutta4;
 
  function f(x, y: real): real;
  begin
    f := y + cos(x) - sin(x);
  end;
 
  function Yetalon(x: real): real;
  begin
    Yetalon := exp(x) + sin(x);
  end;
 
  procedure RungeKutta4Step(var t, y: real; h: real);
  var
    k1, k2, k3, k4: real;
  begin
    k1 := f(t, y);
    k2 := f(t + (h / 2), y + k1 * (h / 2));
    k3 := f(t + (h / 2), y + k2 * (h / 2));
    k4 := f(t + h, y + k3 * h);
    y  := y + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6;
    t  := t + h;
  end;
 
  procedure RungeKutta4(t0, tfin, h, y0: real; var yfin: real; Nprint: integer);
  var
    t: real;
    y: real;
    Xprint: real;
  begin
    t := t0;
    y := y0;
    if Nprint > 0 then
    begin
      Xprint := t + (tfin - t0) / Nprint;
      writeln('x': 9, 'y': 12, 'Yetalon': 15);
      writeln(t: 12: 5, y: 12: 5, Yetalon(t): 12: 5);
    end;
    while t <= tfin do
    begin
      RungeKutta4Step(t, y, h);
      if (Nprint > 0) and (t >= Xprint) then
      begin
        writeln(t: 12: 5, y: 12: 5, Yetalon(t): 12: 5);
        Xprint := Xprint + (tfin - t0) / Nprint;
      end;
    end;
    yfin := y;
  end;
 
  procedure Differ(t0, tfin, y0: real; eps: real; Nprint: integer);
  var
    h: real;
    y1, y2: real;
  begin
    if Nprint <= 0 then
      exit;
    h := (tfin - t0) / Nprint;
    RungeKutta4(t0, tfin, h, y0, y2, 0);
    repeat
      y1 := y2;
      h  := h / 2;
      RungeKutta4(t0, tfin, h, y0, y2, 0);
    until abs(y2 - y1) < eps;
    RungeKutta4(t0, tfin, h, y0, y2, Nprint);
  end;
 
const
  a  = 0;
  b  = 1;
  Ya = 1;
  Eps = 0.1;
  Nprint = 10;
begin
  Differ(a, b, Ya, Eps, Nprint);
end.