{ program ode1
  resi rovnici y'' + y = 0 s pocatecnimi podminkami y(0)=1, y'=(0) = 0
  od 0 do xmax s presnosti eps a reseni porovnava s y = cos(x), y'=-sin(x)}
CONST
   nvar = 10;
   nstepp = 200;
TYPE
   RealArrayNVAR = ARRAY [1..nvar] OF real;
VAR
   OdeintKmax,OdeintKount: integer;
   OdeintDxsav: real;
   OdeintXp: ARRAY [1..nstepp] OF real;
   OdeintYp: ARRAY [1..nvar,1..nstepp] OF real;
   hdid,hnext:real;
   PocVypDeriv : Integer;
PROCEDURE derivs(x: real;
        VAR y,dydx: RealArrayNVAR);
BEGIN
    DYDX[1] := Y[2];
    DYDX[2] := - Y[1];
    PocVypDeriv := PocVypDeriv +1;
END;
PROCEDURE rk4(VAR y,dydx: RealArrayNVAR;
                       n: integer;
                     x,h: real;
                VAR yout: RealArrayNVAR);
VAR
   i: integer;
   xh,hh,h6: real;
   dym,dyt,yt: ^RealArrayNVAR;
BEGIN
   new(dym);
   new(dyt);
   new(yt);
   hh := h*0.5;
   h6 := h/6.0;
   xh := x+hh;
   FOR i := 1 TO n DO
      yt^[i] := y[i]+hh*dydx[i];
   derivs(xh,yt^,dyt^);
   FOR i := 1 TO n DO
      yt^[i] := y[i]+hh*dyt^[i];
   derivs(xh,yt^,dym^);
   FOR i := 1 TO n DO BEGIN
      yt^[i] := y[i]+h*dym^[i];
      dym^[i] := dyt^[i]+dym^[i]
   END;
   derivs(x+h,yt^,dyt^);
   FOR i := 1 TO n DO
      yout[i] := y[i]+h6*(dydx[i]+dyt^[i]+2.0*dym^[i]);
   dispose(yt);
   dispose(dyt);
   dispose(dym)
END;
PROCEDURE rkqc(VAR y,dydx: RealArrayNVAR;
                        n: integer;
                    VAR x: real;
                 htry,eps: real;
                VAR yscal: RealArrayNVAR;
           VAR hdid,hnext: real);
LABEL 99;
CONST
   pgrow = -0.20;
   pshrnk = -0.25;
   fcor = 0.06666666;
   safety = 0.9;
   errcon = 6.0e-4;
VAR
   i: integer;
   xsav,hh,h,temp,errmax: real;
   dysav,ysav,ytemp: ^RealArrayNVAR;
BEGIN
   new(dysav);
   new(ysav);
   new(ytemp);
   xsav := x;
   FOR i := 1 TO n DO BEGIN
      ysav^[i] := y[i];
      dysav^[i] := dydx[i]
   END;
   h := htry;
   WHILE true DO BEGIN
      hh := 0.5*h;
      rk4(ysav^,dysav^,n,xsav,hh,ytemp^);
      x := xsav+hh;
      derivs(x,ytemp^,dydx);
      rk4(ytemp^,dydx,n,x,hh,y);
      x := xsav+h;
      IF x = xsav THEN BEGIN
         writeln('pause in routine RKQC');
         writeln('stepsize too small');
         readln
      END;
      rk4(ysav^,dysav^,n,xsav,h,ytemp^);
      errmax := 0.0;
      FOR i := 1 TO n DO BEGIN
         ytemp^[i] := y[i]-ytemp^[i];
         temp := abs(ytemp^[i]/yscal[i]);
         IF errmax < temp THEN errmax := temp
      END;
      errmax := errmax/eps;
      IF errmax <= 1.0 THEN BEGIN
         hdid := h;
         IF errmax > errcon THEN
            hnext := safety*h*exp(pgrow*ln(errmax))
         ELSE
            hnext := 4.0*h;
         GOTO 99
      END
      ELSE
         h := safety*h*exp(pshrnk*ln(errmax));
   END;
99:
   FOR i := 1 TO n DO
      y[i] := y[i]+ytemp^[i]*fcor;
   dispose(ytemp);
   dispose(ysav);
   dispose(dysav)
END;
PROCEDURE odeint(VAR ystart: RealArrayNVAR;
                          n: integer;
          x1,x2,eps,h1,hmin: real;
               VAR nok,nbad: integer);
LABEL 99;
CONST
   maxstp = 10000;
   tiny = 1.0e-30;
VAR
   nstp,i: integer;
   xsav,x,{hnext,hdid,}h: real;
   yscal,y,dydx: ^RealArrayNVAR;
BEGIN
   new(yscal);
   new(y);
   new(dydx);
   x := x1;
   IF x2 >= x1 THEN h := abs(h1) ELSE h := -abs(h1);
   nok := 0;
   nbad := 0;
   OdeintKount := 0;
   FOR i := 1 TO n DO y^[i] := ystart[i];
   IF OdeintKmax > 0 THEN
      xsav := x-2.0*OdeintDxsav;
   FOR nstp := 1 TO maxstp DO BEGIN
      derivs(x,y^,dydx^);
      FOR i := 1 TO n DO
         yscal^[i] := abs(y^[i])+abs(dydx^[i]*h)+tiny;
      IF OdeintKmax > 0 THEN
         IF abs(x-xsav) > abs(OdeintDxsav) THEN
            IF OdeintKount < OdeintKmax-1 THEN BEGIN
               OdeintKount := OdeintKount+1;
               OdeintXp[OdeintKount] := x;
               FOR i := 1 TO n DO
                  OdeintYp[i,OdeintKount] := y^[i];
               xsav := x
            END;
      IF (x+h-x2)*(x+h-x1) > 0.0 THEN
         h := x2-x;
      rkqc(y^,dydx^,n,x,h,eps,yscal^,hdid,hnext);
      IF hdid = h THEN
         nok := nok+1
      ELSE
         nbad := nbad+1;
{      WRITEln('od',hdid,' ',nok,' ',nbad); Readln;}
      IF (x-x2)*(x2-x1) >= 0.0 THEN BEGIN
         FOR i := 1 TO n DO ystart[i] := y^[i];
         IF OdeintKmax <> 0 THEN BEGIN
            OdeintKount := OdeintKount+1;
            OdeintXp[OdeintKount] := x;
            FOR i := 1 TO n DO
               OdeintYp[i,OdeintKount] := y^[i]
         END;
         GOTO 99
      END;
      IF abs(hnext) < hmin THEN BEGIN
         writeln('pause in routine ODEINT');
         writeln('stepsize too small');
         readln
      END;
      h := hnext;
   END;
   writeln('pause in routine ODEINT - too many steps');
   readln;
99:
   dispose(dydx);
   dispose(y);
   dispose(yscal)
END;
{ Zacatek hlavniho programu }
VAR
    yhod,ystar: RealArrayNVAR;
    xs,xn,eps,h1,hmin,xmax,dx: real;
    nd, nok, nbad, ncyk, i: integer;
    f1,f2 : text;
BEGIN
   OdeintKmax := 200;
   ASSIGN(f2,'ode1.dat');
   REwrite(f2);
   PocVypDeriv := 0;
   nd := 2;
   eps := 1.e-4;
   hmin := 1.e-10;
   h1 := 1.e-4;
   xn := 0.;
   yhod[1] := 1.;
   yhod[2] := 0.;
   WRIteln(xn,yhod[1],yhod[2]);
   writeln(' zadej xmax, ncyk, eps');
   readln( xmax,ncyk,eps) ;
   dx := xmax / ncyk;
   for i := 1 TO ncyk  DO BEGIN
         xs := xn;
         ystar[1] :=yhod[1];
         ystar[2] :=yhod[2];
         xn := xs + dx;
         OdeintDXsav := dx/100;
         odeint(yhod, nd, xs,xn, eps,h1,hmin, nok,nbad);
         writeln(nok:5,nbad:5,' h1 ',h1:10,' hdid ',hdid:10,
                 ' hnext ',hnext:10 );
         WRIteln(f2,'  ',xn:12,'  ',yhod[1]:12,'  ',yhod[2]:12,
                    '  ',cos(xn):12,'  ',-sin(xn):12);
         WRIteln('  ',xn:12,'  ',yhod[1]:12,'  ',yhod[2]:12,
                    '  ',cos(xn):12,'  ',-sin(xn):12);
         if(OdeintKount > 1) then
             h1 := OdeintXP[OdeintKount] - OdeintXP[OdeintKount-1]
         else
             h1 := xn - xs;
         IF (h1 > 0.7*dx) THEN h1 := dx;
         IF i MOD 10=0 THEN readln;
      END;
      Close(f2);
      Writeln('PocVypDeriv = ',PocVypDeriv);
      Writeln('Konec ... Stiskni ENTER');Readln;
END.