[Delphi] Решение ДУ методом Рунге-Кутта

const
  eps = 0.00001;
  x0  = 3;
  y0  = 3;
  type
      TFunc = function(x,y: extended): extended;

function func(x,y: extended): extended; far;
begin
  func := y*cos(x) - 2*sin(2*x);
end;

function RungeKutt(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i   : word;      { счетчик цикла }
  x,h : extended;  { текущая точка и длина шага }
  res : extended;  { переменная для накопления конечного результата функции }
  k1,k2,k3,k4: extended; { вспомогательные переменные вычисления результата }
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  res:= y0_; { устанавливаем начальные значения}
  x:=x0_;
  for i:=1 to n do
    begin { вычисляем результат по методу Рунге-Кутта 4го порядка }
      k1:=f(x,res);
      k2:=f(x+h/2,res+h*k1/2);
      k3:=f(x+h/2,res+h*k2/2);
      k4:=f(x+h,res+h*k3);
      res:=res+h*(k1+2*k2+2*k3+k4)/6;
      x:=x+h; { переходим к следующей точке }
    end;
  RungeKutt:=res; { присваиваем конечный результат функции }
end;