Решение уравнения методом Рунге-Кутты 4-го порядка

Const
 h=0.1;
 N=10;

Var
 X,Y:Array[0..N] of Real;
 K1,K2,K3,K4,Sum:Real;
 i:Word;

Function F(Xf,Yf:Real):Real;
begin
 F:=Cos(Xf*4+Yf)+(Xf-Yf)*7;
end;

Function fK1(Xk,Yk:Real):Real;
begin
 fK1:=F(Xk,Yk);
end;

Function fK2(Xk,Yk,Hk,Q:Real):Real;
begin
 fK2:=F(Xk+Hk/2,Yk+Hk/2*Q);
end;

Function fK3(Xk,Yk,Hk,Q:Real):Real;
begin
 fK3:=F(Xk+Hk/2,Yk+Hk/2*Q);
end;

Function fK4(Xk,Yk,Hk,Q:Real):Real;
begin
 fK4:=F(Xk+Hk,Yk+Hk*Q);
end;

BEGIN
 X[0]:=0;
 Y[0]:=0;
 WriteLn(' X           Y');
 Writeln(X[0]:4:2,Y[0]:15:5);
 For i:=1 to N do
  begin
   X[i]:=X[i-1]+h;
   K1:=fK1(X[i-1],Y[i-1]);
   K2:=fK2(X[i-1],Y[i-1],h,K1);
   K3:=fK3(X[i-1],Y[i-1],h,K2);
   K4:=fK4(X[i-1],Y[i-1],h,K3);
   Y[i]:=Y[i-1]+h/6*(K1+2.0*K2+2.0*K3+K4);
   Writeln(X[i]:4:2,Y[i]:15:5);
  end;
 Sum:=(Sqr(Y[0])+Sqr(Y[N]))/3*h;
 for i:=1 to (N div 2) do
  Sum:=Sum+4.0*h/3.0*Sqr(Y[2*i-1]);
 for i:=2 to (N div 2) do
  Sum:=Sum+2.0*h/3.0*Sqr(Y[2*i-2]);
 Writeln('Heat: ', Sum:0:5);
 Readln
End.