{$N+}{$E+}
{ program ode41 - ukazka reseni stiff rovnice
       u' =  998 u + 1998 v ;        u(0) = 1;
       v' = -999 u - 1999 v ;        v(0) = 0;
  reseni je u = 2 exp(-x) - exp(-1000x);
            v = - exp(-x) + exp(-1000x);
  prestoze druhy clen je pro x > 0.1 zcela zanedbatelny,
  diktuje stale malou delku kroku ----> Stiff-stable metody}
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;
   PocVypDeriv: LongInt;
PROCEDURE derivs(x: real;
        VAR y,dydx: RealArrayNVAR);
BEGIN
    DYDX[1] := 998*Y[1] + 1998*Y[2];
    DYDX[2] := -999*Y[1] - 1999*Y[2];
    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-10;
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);}
      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;
   Function mexp(x:real):real;
BEGIN
   if(x > -84) then
      mexp := exp(x)
   else
      mexp := 0.;
END;
{ Zacatek hlavniho programu }
VAR
    yhod,ystar: RealArrayNVAR;
    xs,xn,eps,h1,hmin,xmax,dx,hdid,hu,hv: real;
    nd, nok, nbad, ncyk, i: integer;
    f1,f2 : text;
BEGIN
   OdeintKmax := 200;
   ASSIGN(f2,'ode41.dat');
   REwrite(f2);
   PocVypDeriv := 0;
   nd := 2;
   eps := 1.e-4;
   hmin := 1.e-7;
   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) ;
   Odeintdxsav := xmax / ncyk;
   xs := xn;
   ystar[1] :=yhod[1];
   ystar[2] :=yhod[2];
   xn := xmax;
   odeint(yhod, nd, xs,xn, eps,h1,hmin, nok,nbad);
   Writeln('PocVypDeriv = ',PocVypDeriv);
         hu := 2*exp(-xn) - mexp(-1000*xn);
         hv := -exp(-xn)  + mexp(-1000*xn);
         writeln(nok:5,nbad:5,'  ',h1:10);
         WRIteln(f2,'  ',xn:12,'  ',yhod[1]:12,'  ',yhod[2]:12,
                '  ',hu:12,'  ',hv:12);
         WRIteln('  ',xn:12,'  ',yhod[1]:12,'  ',yhod[2]:12,
                 '  ', hu:12,' ', hv:12 );
    for i:= 1 to OdeintKount DO BEGIN
         xs := OdeintXp[i];
         ystar[1] := OdeintYp[1,i];
         ystar[2] := OdeintYp[2,i];
         hu := 2*exp(-xs) - mexp(-1000*xs);
         hv := -exp(-xs) + mexp(-1000*xs);
         WRIteln(f2,'  ',xs:12,'  ',ystar[1]:12,'  ',ystar[2]:12,
                 '  ', hu:12,' ', hv:12 );
         WRIteln('  ',xs:14,'  ',ystar[1]:14,'  ',ystar[2]:14,
                 '  ', hu:14,' ', hv:14 );
    END;
    Close(f2);
    Writeln('Koncim... Stikni ENTER'); Readln
END.