Метод Эйлера

PROGRAM dif_ur;
USES
   crt;
CONST
   e = 0.005;         {точность решения}
   a = 0;             {начало отрезка}
   b = 5;             {конец отрезка}
   N = 10;            {число отрезков, на которые поделили [ab]}
   h =(b-a)/N;        {шаг}
   por_sys = 2;       {порядок системы}
TYPE
   type1 = array[1..por_sys] of real;
VAR
   h,x: real;          {глобальные переменные}
   y: type1;           {глобальные переменные}
   vspom: type1;
   i,j,scr_y: integer; {вспомогательная переменная}

Function F(x:real;y:type1;key:integer):real;
  begin
     case key of
        1:f:=(-x*y[1]+y[2])/(x*x+1); {1-e уравнение системы}
        2:f:=(-y[1]-x*y[2])/(x*x+1); {2-e уравнение системы}
     end;
  end;

Function y_PP(x:real;y:type1;key:integer):real;
  begin
     case key of
        1: y_PP:=(x*x*y[1]-y[1]-2*x*y[2])/sqr(x*x+1);
[RIGHT]{производная 1–го уравнения системы}[/RIGHT]
        2: y_PP:=(x*x*y[2]-y[2]+2*x*y[1])/sqr(x*x+1);
[RIGHT]{производная 2–го уравнения системы}[/RIGHT]
  end;
end;

   {погрешность метода Эйлера на одном шаге итерации}
 Function Pogrechnost(h:real):real;
    var
       i: integer;
       Max_pr: real;
    begin
       Max_pr:=y_PP(x,y,1);
       for i:=2 to por_sys do
          if y_PP(x,y,i)>Max_pr
             then
                Max_pr:=y_PP(x,y,i);
       pogrechnost:=(h*h*Max_pr)/2;
    end;

{вывод на экран решения системы дифференциальных уравнений}
 Procedure PRINT;
  begin
     gotoxy (3,2);      writeln (x:5:3);
     gotoxy (23,scr_y); writeln (y[1]:9:4);
     gotoxy (43,scr_y); writeln (y[2]:9:4);
     gotoxy (63,scr_y); writeln (pogrechnost(h):9:4);
     scr_y:=scr_y+1;
  end;

 BEGIN
    clrscr;
    Scr_y:=3;
    x:=a;

    y[1]:=1;
    y[2]:=1;               {условие Коши, начальное условие i=0}

    gotoxy (3,2);  writeln ('Шаг');
    gotoxy (23,2); writeln ('Y1');
    gotoxy (43,2); writeln ('Y2');
    gotoxy (63,2); writeln ('Погрешность');
    while x<b do
       begin
          x:=x+h;
          for j:=1 to por_sys do
             begin
                vspom[j]:=F(x,y,j);
                y[j]:=y[j]+vspom[j]*h;
             end;
          PRINT;
       end;
    readkey
 END.