{$N+}{$E+}
{ program BS41 - ukazka (nevhodna) Bulirsch-Storovy metody
  pro reseni ODE - zde 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;
CONST
   RzextrImax  =  11;
   RzextrNmax  =  10;
   RzextrNcol  =  7;
VAR
   RzextrX: ARRAY [1..RzextrImax] OF real;
   RzextrD: ARRAY [1..RzextrNmax,1..RzextrNcol] OF real;
PROCEDURE rzextr(iest: integer;
                 xest: real;
       VAR yest,yz,dy: RealArrayNVAR;
               n,nuse: integer);
VAR
   m1,k,j: integer;
   yy,v,ddy,c,b1,b: real;
   fx: ARRAY [1..RzextrNcol] OF real;
BEGIN
   RzextrX[iest] := xest;
   IF iest = 1 THEN
      FOR j := 1 TO n DO BEGIN
         yz[j] := yest[j];
         RzextrD[j,1] := yest[j];
         dy[j] := yest[j]
      END
   ELSE BEGIN
      IF iest < nuse THEN
         m1 := iest
      ELSE
         m1 := nuse;
      FOR k := 1 TO m1-1 DO
         fx[k+1] := RzextrX[iest-k]/xest;
      FOR j := 1 TO n DO BEGIN
         yy := yest[j];
         v := RzextrD[j,1];
         c := yy;
         RzextrD[j,1] := yy;
         FOR k := 2 TO m1 DO BEGIN
            b1 := fx[k]*v;
            b := b1-c;
            IF b <> 0.0 THEN BEGIN
               b := (c-v)/b;
               ddy := c*b;
               c := b1*b
            END
            ELSE
               ddy := v;
            IF k <> m1 THEN
               v := RzextrD[j,k];
            RzextrD[j,k] := ddy;
            yy := yy+ddy
         END;
         dy[j] := ddy;
         yz[j] := yy
      END
   END
END;
PROCEDURE mmid(VAR y,dydx: RealArrayNVAR;
                        n: integer;
                  xs,htot: real;
                    nstep: integer;
                 VAR yout: RealArrayNVAR);
VAR
   step,i: integer;
   x,swap,h2,h: real;
   ym,yn: ^RealArrayNVAR;
BEGIN
   new(ym);
   new(yn);
   h := htot/nstep;
   FOR i := 1 TO n DO BEGIN
      ym^[i] := y[i];
      yn^[i] := y[i]+h*dydx[i]
   END;
   x := xs+h;
   derivs(x,yn^,yout);
   h2 := 2.0*h;
   FOR step := 2 TO nstep DO BEGIN
      FOR i := 1 TO n DO BEGIN
         swap := ym^[i]+h2*yout[i];
         ym^[i] := yn^[i];
         yn^[i] := swap
      END;
      x := x+h;
      derivs(x,yn^,yout)
   END;
   FOR i := 1 TO n DO
      yout[i] := 0.5*(ym^[i]+yn^[i]+h*yout[i]);
   dispose(yn);
   dispose(ym)
END;

PROCEDURE bsstep(VAR y,dydx: RealArrayNVAR;
                          n: integer;
                      VAR x: real;
                   htry,eps: real;
                  VAR yscal: RealArrayNVAR;
             VAR hdid,hnext: real);
LABEL 99;
CONST
   imax = 11;
   nuse = 7;
   shrink = 0.95e0;
   grow = 1.2e0;
VAR
   j,i: integer;
   xsav,xest,h,errmax: real;
   ysav,dysav,yseq,yerr: ^RealArrayNVAR;
   nseq: ARRAY [1..imax] OF integer;
BEGIN
   new(ysav);
   new(dysav);
   new(yseq);
   new(yerr);
   nseq[1] := 2;
   nseq[2] := 4;
   nseq[3] := 6;
   nseq[4] := 8;
   nseq[5] := 12;
   nseq[6] := 16;
   nseq[7] := 24;
   nseq[8] := 32;
   nseq[9] := 48;
   nseq[10] := 64;
   nseq[11] := 96;
   h := htry;
   xsav := x;
   FOR i := 1 TO n DO BEGIN
      ysav^[i] := y[i];
      dysav^[i] := dydx[i]
   END;
   WHILE true DO BEGIN
      FOR i := 1 TO imax DO BEGIN
         mmid(ysav^,dysav^,n,xsav,h,nseq[i],yseq^);
         xest := sqr(h/nseq[i]);
         rzextr(i,xest,yseq^,y,yerr^,n,nuse);
         IF i > 3 THEN BEGIN
            errmax := 0.0;
            FOR j := 1 TO n DO
               IF errmax < abs(yerr^[j]/yscal[j]) THEN
                  errmax := abs(yerr^[j]/yscal[j]);
            errmax := errmax/eps;
            IF errmax < 1.0 THEN BEGIN
               x := x+h;
               hdid := h;
               IF i = nuse THEN
                  hnext := h*shrink
               ELSE IF i = nuse-1 THEN
                  hnext := h*grow
               ELSE
                  hnext := (h*nseq[nuse-1])/nseq[i];
               GOTO 99
            END
         END
      END;
      h := 0.25*h;
      FOR i := 1 TO (imax-nuse) DIV 2 DO h := h/2;
      IF x+h = x THEN BEGIN
         writeln('pause in routine BSSTEP');
         writeln('step size underflow');
         readln
      END
   END;
99:
   dispose(yerr);
   dispose(yseq);
   dispose(dysav);
   dispose(ysav)
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;
      BSSTEP(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,'bsode41.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:12,'  ',ystar[1]:12,'  ',ystar[2]:12,
                 '  ', hu:12,' ', hv:12 );
    END;
    Close(f2);
    Writeln('Konec ... Stiskni ENTER');Readln;
END.