Program Miln;

Uses Crt;

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

  Clrscr;

  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.