[Страдалецъ]

Код:

program Miln;

{$APPTYPE CONSOLE}

uses
  SysUtils;

var
  x0,y0,z0,h:real; 
  k1,k2,k3,k4:real; 
  r1,r2,r3,r4:real; 
  eps,abs_pogr:real; 
  
  z,zkor,zpr,ypr,ykor,x,y:array [0..10] of real; 
  i:integer; 
  
function f1(xa,ya,za:real):real; 
begin 
  f1:=2*sqr(xa)+2*ya+za; 
end; 
  
function f2(xa,ya,za:real):real; 
begin 
  f2:=1-2*sqr(xa)+2*ya-za; 
end; 

begin
  eps:=10e-6;
  x0:=0; 
  y0:=1; 
  z0:=1; 
  h:=0.1; 
  x[0]:=x0; 
  y[0]:=y0; 
  z[0]:=y0; 
  i:=0; 
  { Gotovim 1-e 3 tochki po metodu runge-kutta } 
  while i<=3  do
    begin 
      k1:=h*f1(x[i],y[i],z[i]); 
      r1:=h*f2(x[i],y[i],z[i]); 
      k2:=h*f1(x[i]+h/2,y[i]+k1/2,z[i]+r1/2); 
      r2:=h*f2(x[i]+h/2,y[i]+k1/2,z[i]+r1/2); 
      k3:=h*f1(x[i]+h/2,y[i]+k2/2,z[i]+r2/2); 
      r3:=h*f2(x[i]+h/2,y[i]+k2/2,z[i]+r2/2); 
      k4:=h*f1(x[i]+h,y[i]+k3,z[i]+r3); 
      r4:=h*f2(x[i]+h,y[i]+k3,z[i]+r3); 
  
      y[i+1]:=y[i]+(k1+2*k2+2*k3+k4)/6; 
      z[i+1]:=z[i]+(r1+2*r2+2*r3+r4)/6; 
  
      x[i+1]:=x[i]+h; 
      i:=i+1;
    end; 
  
  i:=4; 
  while x[i]<=1.0+h{eps> abs_pogr} do
    begin 
  { etap prognoza i korrektsii} 
      ypr[i]:=y[i-4]+(4*h)/3*(2*f1(x[i-3],y[i-3],z[i-3])-f1(x[i-2],y[i-2],z[i-2])+2*f1(x[i-1],y[i-1],z[i-1])); 
  
      ykor[i]:=y[i-2]+(h/3)*(f1(x[i-2],y[i-2],z[i-2])+4*f1(x[i-1],y[i-1],z[i-1])+f1(x[i],ypr[i],z[i])); 

      zpr[i]:=z[i-4]+(4*h)/3*(2*f2(x[i-3],y[i-3],z[i-3])-f2(x[i-2],y[i-2],z[i-2])+2*f2(x[i-1],y[i-1],z[i-1])); 
  
      zkor[i]:=z[i-2]+(h/3)*(f2(x[i-2],y[i-2],z[i-2])+4*f2(x[i-1],y[i-1],z[i-1])+f2(x[i],zpr[i],z[i])); 
  
  
      abs_pogr:=abs(ykor[i]-ypr[i])/29; 
      if abs_pogr>eps then 
        y[i]:=ykor[i] 
      else 
          y[i]:=ypr[i]; 
      abs_pogr:=abs(zkor[i]-zpr[i])/29; 
      if abs_pogr>eps then 
        z[i]:=zkor[i] 
      else 
          z[i]:=zpr[i];
      x[i+1]:=x[i]+h; 
      i:=i+1; 
    end; 
  WriteLn; 
  for i:=0 to 10 do 
    begin 
      WriteLn(x[i]:10:4,' ',y[i]:10:4,' ',z[i]:10:4); 
    end; 
  Readln; 
end.