Решение дифференциального уравнения методом Рунге-Кутта-Мерсона

uses crt,graph;
var x0,y0,gm,gd,i,n:integer;
mx,my,s,h,a,b:real;
x,y:array[0..100] of real;
function f(y,x:real):real;
begin
f:=x+y;
end;
begin
writeln('Vvedit interval [a;b]');
readln(a,b);
writeln('Vvedit po4atkovi umovy');
readln(y[0]);
x[0]:=a;
writeln('Vvedit krok h');
readln(h);
n:=round((b-a)/h);
for i:=1 to n do
begin
s:=h*f(y[i-1],x[i-1]);
y[i]:=s;
s:=h*f(y[i-1]+s/2,x[i-1]+h/2);
y[i]:=y[i]+2*s;
s:=h*f(y[i-1]+s/2,x[i-1]+h/2);
y[i]:=y[i]+2*s;
s:=h*f(y[i-1]+s,x[i-1]+h/2);
y[i]:=y[i]+s;
y[i]:=y[i-1]+y[i]/6;
x[i]:=x[i-1]+h;
end;
for i:=0 to n do
writeln(x[i]:2:3,' ',y[i]:2:3);
readln;
gd:=detect;
initgraph(gd,gm,'');
if graphresult<>grok then halt(1);
a:=y[0];
b:=y[0];
for i:=1 to n do
begin
if a<y[i] then a:=y[i];
if b>y[i] then b:=y[i];
end;
mx:=(getmaxx-20)/(x[n]-x[0]);
my:=(getmaxy-20)/(a-b);
y0:=round(a*my+10);
x0:=round(-x[0]*mx+10);
line(x0,0,x0,getmaxy);
line(0,y0,getmaxx,y0);
setcolor(red);
for i:=1 to n do
line(round(x[i]*mx-x0),round(-y[i]*my+y0),round(x[i-1]*mx-x0),round(-y[i-1]*my+y0));
readln;
closegraph;
writeln(b);
end.