Решение дифференциального уравнения методом Рунге-Кутта четвертого порядка - Pascal

program RungeKutta4;
 
  function f(x, y: real): real;
  begin
    f := x * sin(x) + y / x;
  end;
 
  function Yetalon(x: real): real;
  begin
    Yetalon := 2 * x - x * cos(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;
    Tprint: real;
  begin
    t := t0;
    y := y0;
    if Nprint > 0 then
    begin
      Tprint := 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 + h) <= (tfin + h / 2) do
    begin
      RungeKutta4Step(t, y, h);
      if (Nprint > 0) and (abs(t - Tprint) < h / 2) then
      begin
        writeln(t: 12: 5, y: 12: 5, Yetalon(t): 12: 5);
        Tprint := Tprint + (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  = Pi / 2;
  b  = Pi / 2 + 1;
  Ya = Pi;
  Eps = 0.01;
  Nprint = 10;
var
  yfin: real;
begin
  writeln('Решение с фиксированным шагом ', (b - a) / Nprint: 0: 5);
  RungeKutta4(a, b, (b - a) / Nprint, Ya, yfin, Nprint);
  writeln('Решение с заданной точностью Eps=', Eps: 0: 5);
  Differ(a, b, Ya, Eps, Nprint);
end.