Метод Адамса 4-го порядка

PROGRAM Ad_2;
  USES
     CRT;
  LABEL
     MET;
  CONST
     a = 1;
     b = 2;
     e = 0.01;
  VAR
     h : Real;
     i,u : Integer;
     Y1,Y2,Y3,X : Array [1..30] of Real;
  BEGIN
    clrscr;
    h:=0.05;
    Y1[1]:=1;
    Y2[1]:=1;
    Y3[1]:=1;
    repeat
          u:=0;
MET:   for i:=1 to 2 do
          begin
             Y1[i+1]:=Y1[i]+h*(Y1[i]+h/2*(3*Y2[i]-2*Y3[i]));
             Y2[i+1]:=Y2[i]+h*(Y2[i]+h/2*(-Y3[i]-3*Y1[i]));
             Y3[i+1]:=Y3[i]+h*(Y3[i]+h/2*(2*Y1[i]-Y2[i]));
          end;
       if u=0
          then
             begin
                h:=h/2;
                u:=1;
                goto MET
             end;
     until ((abs(Y1[2])-Y1[3])<=(3*e)) and ((abs(Y2[2])-Y2[3])<=(3*e))
                                       and ((abs(Y3[2])-Y3[3])<=(3*e));
    h:=h*2;
    for i:=2 to (round((b-a)/h)+1) do
       begin
          Y1[i+1]:=Y1[i]+h/2*(3*(3*Y2[i]-2*Y3[i])-(3*Y2[i-1]-2*Y3[i-1]));
          Y2[i+1]:=Y2[i]+h/2*(3*(-Y3[i]-3*Y1[i])-(-Y3[i-1]-3*Y1[i-1]));
          Y3[i+1]:=Y3[i]+h/2*(3*(2*Y1[i]-Y2[i])-(2*Y1[i-1]-Y2[i-1]))
       end;
    writeln ('   X           Y1          Y2          Y3');
    writeln ('-------------------------------------------');
    for i:=1 to (round((b-a)/h)+1) do
       begin
          X[i]:=a+h*(i-1);
          writeln (X[i]:7:4, Y1[i]:12:4, Y2[i]:12:4, Y3[i]:12:4);
       end;
    writeln;
    write (' ':5,'Шаг интегрирования: ',h:4:4);
    readkey
  END.

program Project1; 
 
uses 
 crt; 
 
var 
 x,y,x0,xn,y0,h, 
 k1,k2,k3,k4:real; 
 n,i: integer; 
 Fi : array [1..30]of real;
 
function fxy(x,y: real): real; 
begin 
 fxy:=x+y;
end; 
 
begin 
 Write ('Vvedite diapazon otrezka (a,b): '); 
 Readln(x0,xn); 
 Write ('Vvedite y(', x0:2:2, '): '); 
 Readln (y0); 
 Write ('Vvedite kolichestvo shagov izmenenia argumenta: '); 
 Readln (n); 
 h := (xn - x0) / n; 
 x := x0; 
 y := y0;
 {SetLength(Fi[i],n+1);}
 Fi[i] := fxy(x,y);
 
 Writeln('Nachalnoe znachenie: ', ' x = ', x:4:4, ', y = ', y:4:4); 
 
 {//Нахождение первых 3-ёх точек (4-тая f0) методом Рунге-Кутта }
 for i := 1 to 3 do 
 begin 
 k1 := h*fxy(x,y); 
 Write(k1:4:4); 
 k2 := h*fxy(x+h/2,y+k1/2); 
 k3 := h*fxy(x+h/2,y+k2/2); 
 k4 := h*fxy(x+h,y+k3); 
 y := y + h * (k1+2*k2+2*k3+k4) / 6; 
 x := x + h; 
 Fi[i]:= fxy(x,y);
 Writeln('Shag #', i, ' , x = ', x:4:4, ', y = ', y:4:4); 
 readkey; 
 end; 
 
 {//Нахождение остальных точек методом Адамса  }
 for i := 4 to n do 
 begin 
 y := y + h * (55*Fi[i-1]-59*Fi[i-2]+37*Fi[i-3]-9*Fi[i-4])/24; 
 x := x + h; 
 Fi [i]:= fxy(x,y);
 Writeln('Shag #', i, ' , x = ', x:4:4, ', y = ', y:4:4); 
 readkey; 
 end; 
 
{//Проверка другим методом }
 x := x0; 
 y := y0; 
 for i := 1 to n do 
 begin 
 k1 := h*fxy(x,y); 
 k2 := h*fxy(x+h/2,y+k1/2); 
 k3 := h*fxy(x+h/2,y+k2/2); 
 k4 := h*fxy(x+h,y+k3); 
 y := y + h * (k1+2*k2+2*k3+k4) / 6; 
 x := x + h; 
 Writeln('Shag #', i, ' , x = ', x:4:4, ', y = ', y:4:4); 
 readkey; 
 end; 
end.