Thursday, October 12, 2017

Stormer's method

'*************************************************************
'*   Differential equation y"=f(x,y,y') by Stormer's method  *
'* --------------------------------------------------------- *
'* SAMPLE RUN:                                               *
'* Integrate y" = 8yy / (1 + 2x) from x=0 to x=1,            *
'* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2      *
'* and compare with exact solution: y = 1 / (1 + 2x)         *
'*                                                           *
'* Output file (stormer.lst):                                *
'*                                                           *
'* --------------------------------------------------------- *
'*   Differential equation y"=f(x,y,y') by Stormer's method  *
'* --------------------------------------------------------- *
'*        X           Y           Y exact         Error      *
'*     0.000       1.000000       1.000000     0.0000000000  *
'*     0.010       0.980392       0.980392     0.0000000001  *
'*     0.020       0.961538       0.961538     0.0000000295  *
'*     0.030       0.943396       0.943396     0.0000000457  *
'*     0.040       0.925926       0.925926     0.0000000974  *
'*     0.050       0.909091       0.909091     0.0000001285  *
'*     ...         ...            ...          ...           *
'*     0.950       0.344866       0.344828     0.0000381695  *
'*     0.960       0.342505       0.342466     0.0000388874  *
'*     0.970       0.340176       0.340136     0.0000396196  *
'*     0.980       0.337878       0.337838     0.0000403406  *
'*     0.990       0.335612       0.335570     0.0000410721  *
'*     1.000       0.333375       0.333333     0.0000418231  *
'*                                                           *
'*  End of file.                                             *
'*                              Basic Version By J-P Moreau  *
'*                                   (www.jpmoreau.fr)       *
'*************************************************************
defint i-n
defdbl a-h,o-z

dim c(4),x(4),y(4),z(4)

  h=0.01#
  a1=1.08333333333333#
  a2=-2#*(a1-1#)
  a3=a1-1#

  f$="  ##.####   ##.######   ##.######   ##.##########"

  open "stormer.lst" for output as #2

  cls
  print #2, "-----------------------------------------------------------------"
  print #2, "  Differential equation y""=f(x,y,y') by Stormer's method"
  print #2, "-----------------------------------------------------------------"
  'initial conditions
  x(1)=0# : y(1)=1# : z(1)=-2#
  xx=x(1) : gosub 1200 : yex=Fx : er=0#
  print #2,"      X         Y        Y exact         Error"
  print #2, using f$; x(1); y(1); yex; er
  for k=1 to 2
    'call Runge-Kutta for first 2 steps
    gosub 2000
    xx=x(k+1) : gosub 1200 : yex=Fx : er=abs(yex-y(k+1))
    print #2, using f$; x(k+1); y(k+1); yex; er
  next k
'main Stormer loop
10 'continue
    for k=2 to 4
      xx=x(5-k) : yy=y(5-k) : gosub 1100
      c(k)=G
    next k
    y(4)=2#*y(3)-y(2)+h*h*(a1*c(2)+a2*c(3)+a3*c(4))
    x(4)=x(3)+h : xx=x(4) : gosub 1200 : yex=Fx : er=abs(yex-y(4))
    print #2, using f$; x(4); y(4); yex; er
    for k=1 to 3
      x(k)=x(k+1) : y(k)=y(k+1)
    next k
  if x(3)<1# then goto 10  'end x value = 1
  print #2,
  print #2, " End of file."
  Close #2
  print
  print " Results in file stormer.lst"
  print

END

1000 'real*8 function F(x,y,z)
  F=zz
return

1100 'real*8 function G(x,y,z)
  G=8#*yy*yy/(1#+2#*xx)
return

1200 'exact solution real*8 function Fx(x)
  Fx=1#/(1#+2#*xx)
return

2000 'SUBROUTINE RK4D2(x,y,z,h,x1,y1,z1)
  xx=x(k) : yy=y(k) : zz=z(k)
  gosub 1000 : c1=F
  gosub 1100 : d1=G
  xx=x(k)+h/2# : yy=y(k)+h/2#*c1 : zz=z(k)+h/2#*d1
  gosub 1000 : c2=F
  gosub 1100 : d2=G
  zz=z(k)+h/2#*d2 : gosub 1000 : c3=F
  yy=y(k)+h/2#*c2 : gosub 1100 : d3=G
  xx=x(k)+h
  zz=z(k)+h*d3 : gosub 1000 : c4=F
  yy=y(k)+h*c3 : gosub 1100 : d4=G
  x(k+1)=x(k)+h
  y(k+1)=y(k)+h*(c1+2#*c2+2#*c3+c4)/6#
  z(k+1)=z(k)+h*(d1+2#*d2+2#*d3+d4)/6#
return

'End of file stormer.bas
            EXPLANATION FILE OF PROGRAM STORMER
            ===================================



  Solve Y"=f(x,y,y') with initial conditions by Stormer's method
  --------------------------------------------------------------

    The differential equation of order 2 can be replaced by a system of two
  equations of order 1:

  Given y"=f(x,y,y') with y(a) and y'(a), by calling u = y', the problem becomes

                    | u'=f(x,y,u)
                    | y'=u
                    | with y(a), u(a) given


    We start from a particular form of the Taylor's formula, where the remainder
  is under the form of an integral:

                                     x+h
           y(x+h) = y(x) + h y'(x) + Sum (x+h-t) y"(t) dt
                                      x

  In the same way

                                     x-h
           y(x-h) = y(x) - h y'(x) + Sum (x-h-t) y"(t) dt
                                      x

  By summing

                                     x+h       x-h
           y(x+h) - 2y(x) + y(x-h) = Sum ... + Sum ...
                                      x         x

  For the second integral, we change the variable: u = 2x - t

  So       x-h                      x+h
           Sum (x-h-t) y"(t) dt = - Sum (x+h-u) y"(u) du
            x                        x

                                    x+h
           du=-dt ==>           =   Sum (x+h-t) y"(2x-t) dt
                                     x

  Finally
                                     x+h
           y(x+h) - 2y(x) + y(x-h) = Sum (x+h-t)[y"(t)+y"(2x-t)] dt
                                      x

  We now use the interpolation polynomial recalling that y"(x)=f(x,y,y'):

                              x
                               n+1
           y    - y  + y    = Sum  (x   - t)[P(t)+P(2x - 1)] dt
            n+1    n    n-1   x      n+1              n
                               n
              x
            2  n+1                 0             1            2
         = h  Sum (x   - t) [O0 Div (fn) + O1 Div (fn) +O2 Div (fn)] dt
              x     n+1
               n
                       m  1        |(-s)   (m)|
         with O  = (-1)  Sum (1-s) |(  ) + ( )| ds
               m          0        |(m )   (s)|

                                 (see file Adambash.txt).  


  This leads to Stormer's formulas:


    Explicit:  y    - 2 y  + y    = (h²/12)[13f - 2 f    + f   ]
                n+1      n    n-1              n     n-1    n-2

    Implicit:  y    - 2 y  + y    = (h²/12)[f   + 10 f  + f   ]
                n+1      n    n-1            n+1      n    n-1


  From [BIBLI 04]
---------------------------------------
End of file Stormer.txt

Tuesday, October 10, 2017

Implicit Gear method of order 4

DECLARE SUB SWAP1 (a#, b#)
DECLARE FUNCTION det# (n%, xmat#())
DECLARE FUNCTION XNormMax# (vektor#(), n%)
DECLARE SUB settxt (bspn AS INTEGER, n%, dgltxt$())
DECLARE SUB gear4 (x#, xend#, bspn AS INTEGER, n%, y#(), epsabs#, epsrel#, h#, 
fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER)
DECLARE FUNCTION DistMax# (vector1#(), vector2#(), n%)
DECLARE SUB ruku23 (x#, y#(), bspn AS INTEGER, n%, h#, y2#(), y3#())
DECLARE SUB dgl (bspn AS INTEGER, n%, x#, y#(), f#())
DECLARE SUB engl45 (x#, y#(), bspn AS INTEGER, n%, h#, y4#(), y5#())
DECLARE SUB prdo45 (x#, y#(), bspn AS INTEGER, n%, h#, y4#(), y5#(), 
steif1 AS INTEGER, steifanz AS INTEGER, steif2 AS INTEGER)
DECLARE SUB awp (x#, xend#, bspn AS INTEGER, n%, y#(), epsabs#, epsrel#, h#, 
methode%, fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER)
DECLARE FUNCTION XMax# (a#, b#)
DECLARE FUNCTION XMin# (a#, b#)
DECLARE SUB ISWAP (ia%, ib%)
DECLARE SUB gauss (mode%, n%, xmat#(), xlumat#(), iperm%(), b#(), x#(), 
signd AS INTEGER, code AS INTEGER)
DECLARE SUB gaudec (n%, xmat#(), xlumat#(), iperm%(), signd AS INTEGER, rc AS INTEGER)
DECLARE SUB gausol (n%, xlumat#(), iperm%(), b#(), x#(), rc AS INTEGER)
'************************************************************************
'*   Solve a first order Stiff System of Differential Equations using   *
'*   the implicit Gear method of order 4.                               *
'* -------------------------------------------------------------------- *
'* Mode of operation:                                                   *
'* ==================                                                   *
'* This program can solve two of the 11 examples of file t_dgls using   *
'* the implicit Gear method of order 4 (see file gear.pas).             *
'* To test other systems of DEs, please proceed as explained in file    *
'* t_dgls.pas.                                                          *
'*                                                                      *
'*   Inputs:                                                            *
'*   =======                                                            *
'*   bspnummer  Number of DE system from t_dgls.pas                     *
'*   epsabs     desired absolute error bound                            *
'*   epsrel     desired relative error bound                            *
'*   x0         left edge of integration                                *
'*   y0[0]   \  known approximation for the solution at x0              *
'* ..  .      >                                                         *
'*   y0[n-1] /                                                          *
'*   h          initial step size                                       *
'*   xend       right endpoint of integration                           *
'*   fmax       maximal number of calls of the right hand side          *
'*                                                                      *
'*   The size n of the DE system is passed on from t_dgls.bas.          *
'* -------------------------------------------------------------------- *
'* SAMPLE RUN                                                           *
'*                                                                      *
'* Example #1:                                                          *
'* (Solve set of differential equations (n=2):                          *
'*     f[0] = y[0] * y[1] + COS(x) - HALF * SIN(TWO * x);               *
'*     f[1] = y[0] * y[0] + y[1] * y[1] - (ONE + SIN(x));               *
'*  Find values of f(0), f(1) at x=1.5).                                *
'*                                                                      *
'* Input example number (0 to 11): 0                                    *
'* abs. epsilon: 1e-6                                                   *
'* rel. epsilon: 1e-8                                                   *
'* x0: 0                                                                *
'* y0[0]: 0.5                                                           *
'* y0[1]: 0.5                                                           *
'* initial step size h: 0.0001                                          *
'* right edge xend: 1.5                                                 *
'* maximal number of calls of right hand side: 6000                     *
'*                                                                      *
'* Input data:                                                          *
'* -----------                                                          *
'* Example # 0                                                          *
'* Number of DEs = 2                                                    *
'* x0       =  0                                                        *
'* xend     =  1.5                                                      *
'* epsabs   =  0.000001                                                 *
'* epsrel   =  0.00000001                                               *
'* fmax     =  6000                                                     *
'* h        =  0.0001                                                   *
'* y0(0)    =  0.5                                                      *
'* y0(1)    =  0.5                                                      *
'*                                                                      *
'* Output data:                                                         *
'* ------------                                                         *
'* error code from gear4: 0                                             *
'* final local step size: 6.06783655109067E-02                          *
'* number of calls of right hand side: 360                              *
'* Integration stopped at x = 1.5                                       *
'*                                                                      *
'* approximate solution y1(x) =  1.23598612893281                       *
'* approximate solution y2(x) = -0.104949617987246                      *
'*                                                                      *
'* Example #2:                                                          *
'* (Solve set of differential equations (n=5):                          *
'*   f[0] = y[1];                                                       *
'*   f[1] = y[2];                                                       *
'*   f[2] = y[3];                                                       *
'*   f[3] = y[4];                                                       *
'*   f[4] = ((REAL)45.0 * y[2] * y[3] * y[4] -                          *
'*          (REAL)40.0 * y[3] * y[3] * y[3]) / (NINE * y[2] * y[2]);    *
'*  Find values of f(0), ..., f(4) at x=1.5).                           *
'*                                                                      *
'* Input example number (0 to 11): 3                                    *
'* abs. epsilon: 1e-10                                                  *
'* rel. epsilon: 1e-10                                                  *
'* x0: 0                                                                *
'* y0[0]: 1                                                             *
'* y0[1]: 1                                                             *
'* y0[2]: 1                                                             *
'* y0[3]: 1                                                             *
'* y0[4]: 1                                                             *
'* initial step size h: 0.001                                           *
'* right edge xend: 1.5                                                 *
'* maximal number of calls of right hand side: 6000                     *
'*                                                                      *
'* Input data:                                                          *
'* -----------                                                          *
'* Example # 3                                                          *
'* Number of DEs = 5                                                    *
'* x0       =  0                                                        *
'* xend     =  1.5                                                      *
'* epsabs   =  0.0000000001                                             *
'* epsrel   =  0.0000000001                                             *
'* fmax     =  6000                                                     *
'* h        =  0.001                                                    *
'* y0[0]    =  1                                                        *
'* y0[1]    =  1                                                        *
'* y0[2]    =  1                                                        *
'* y0[3]    =  1                                                        *
'* y0[4]    =  1                                                        *
'*                                                                      *
'* Output data:                                                         *
'* ------------                                                         *
'* error code from gear4: 0                                             *
'* final local step size: 4.86347661993806E-03                          *
'* number of calls of right hand side: 3423                             *
'* Integration stopped at x = 1.5                                       *
'*                                                                      *
'* approximate solution y1(x) =  4.36396102990278                       *
'* approximate solution y2(x) =  4.00000000763431                       *
'* approximate solution y3(x) =  2.82842715661993                       *
'* approximate solution y4(x) =  4.86163232512205E-08                   *
'* approximate solution y5(x) = -3.7712362229557                        *
'*                                                                      *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
'*                                                                      *
'*                            Quick Basic Release By J-P Moreau, Paris. *
'*                                       (www.jpmoreau.fr)              *
'************************************************************************
DefDbl A-H, O-Z
DefInt I-N

Option Base 0
  
' epsabs,                  absolute error bound
' epsrel,                  relative error bound
' x0                       left edge of integration interval
' y0 (0:n-1)-vector        initial values
' yex(0:n-1)-vector        exact solution (when given)
' h                        initial, final step size
' xend                     right edge of integration interval

Dim fmax As Integer      ' maximal number of calls of right side in gear4
Dim aufrufe As Integer   ' actual number of function calls
Dim bspnummer As Integer ' # example
' n                      ' number of DEs in system (see t_dlgs.bas)
Dim fehler As Integer    ' error code from subroutine gear4
Dim dgltxt$(5)           ' text of equations
' -------------------- read input  --------------------
  Cls
  Print
  INPUT " Example number (0 to 11): ", bspnummer
  If bspnummer <> 0 And bspnummer <> 3 Then
    Print
    Print " Example not registered."
    End   'stop program
  End If

' input absolute and relative errors
  INPUT " epsabs = ", epsabs
  INPUT " epsrel = ", epsrel

' input x starting value
  INPUT " x0 = ", x0

  Call settxt(bspnummer, n, dgltxt$())     'read text of equations
                                           'and parameter n
' allocate memory for vectors y0, yex
  Dim y0(n)
  Dim yex(n)

' input initial values y0(i)
' y0(0) = 0.5: y0(1) = 0.5 (Ex. 0)
  For i = 0 To n - 1
    PRINT " y0("; i; ") = "; : INPUT "", y0(i)
  Next i

' input initial step size
  INPUT " h = ", h

' input ending x value
  INPUT " xend = ", xend
' input maximum number of calls to right hand side
  INPUT " fmax = ", fmax

' ----------------- print input data ----------------------
  Cls
  Print " =============================================="
  Print "   Solve a first order ordinary system of DEs  "
  Print " using the implicit method of Gear of 4th order"
  Print " =============================================="
  Print
  Print " System of DEs:"
  Print " ------------- "
  For i = 0 To n - 1
    Print " "; dgltxt$(i)
  Next i
  Print " Example # "; bspnummer
  Print " Number of DEs = "; n
  Print " x0     = "; x0
  Print " xend   = "; xend
  Print " epsabs = "; epsabs
  Print " epsrel = "; epsrel
  Print " fmax   = "; fmax
  Print " h      = "; h
  For i = 0 To n - 1
    Print " y0("; i; ") = "; y0(i)
  Next i

' ------------------- Solve system of DEs ------------------------------
  gear4 x0, xend, bspnummer, n, y0(), epsabs, epsrel, h, fmax, aufrufe, fehler

  If fehler <> 0 Then     'something went wrong
    Print " Gear4: error nø "; 10 + fehler
    End   'stop program
  End If

' ---------------------- print results ---------------------
  Print
  Print " Output data:"
  Print "  -----------"
  Print " error code from gear4: "; fehler
  Print " final local step size: "; h
  Print " number of calls of right hand side: "; aufrufe
  Print " Integration stopped at x = "; x0
  Print
  For i = 0 To n - 1
    Print " approximate solution y"; i + 1; "(x) = "; y0(i)
  Next i
  Print

End 'of main program

'************************************************************************
'*                                                                      *
'* Solve an ordinary system of first order differential equations using *
'* -------------------------------------------------------------------- *
'* automatic step size control                                          *
'* ----------------------------                                         *
'*                                                                      *
'* Programming language: ANSI C                                         *
'* Author:               Klaus Niederdrenk (FORTRAN)                    *
'* Adaptation:           Juergen Dietel, Computer Center, RWTH Aachen   *
'* Source:               existing C, Pascal, QuickBASIC and FORTRAN     *
'*                       codes                                          *
'* Date:                 6.2.1992, 10.2.1995                            *
'*                                                                      *
'*                       Quick Basic Release By J-P Moreau, Paris.      *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996".                      *
'************************************************************************
' 1st order DESs with automatic step size control ...........................
Sub awp(x, xend, bspn As Integer, n, y(), epsabs, epsrel, h, methode, fmax As Integer, aufrufe As Integer, fehler As Integer)
'  double x,                      initial/final x value ..............
'         xend                    desired end point ..................
' integer bspn,                   # example
'         n                       number of DEs ......................
'  double y(0:n-1),               initial/final y value ..............
'         epsabs,                 absolute error bound ...............
'         epsrel,                 relative error bound ...............
'         h                       initial/final step size ............
' integer methode,                desired method (3, 6, 7) ...........
'         fmax,                   maximal # of calls of  dgl() .......
'         aufrufe,                actual # of calls of  dgl() ........
'         fehler                  error code .........................
'************************************************************************
'* Compute the solution y of a system of first order ordinary           *
'* differential equations       y' = f(x,y)   at xend from the given    *
'* initial data (x0, y0).                                               *
'* We use automatic step size control internally so that the error of   *
'* y (absolutely or relatively) lies within the given error bounds      *
'* epsabs and epsrel.                                                   *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x        initial value for x                                         *
'* y        initial values for y(0:n-1)                                 *
'* bspn     # example                                                   *
'* n        number of differential equations                            *
'* dgl      function that evaluates the right hand side of the system   *
'*          y' = f(x,y)  (see t_dgls) (Removed from list of parameters) *
'* xend     end of integration; xend > x0                               *
'* h        initial step size                                           *
'* epsabs   absolute error bound; >= 0; if = 0 we only check the        *
'*          relative error.                                             *
'* epsrel   relative error bound; >= 0; if = 0 we check only the        *
'*          absolute eror.                                              *
'* fmax     max number of evaluations of right hand side in dgl()       *
'* methode  chooses the method                                          *
'*          = 3: Runge-Kutta method of 2nd/3rd order                    *
'*          = 6: England formula of 4th/5th order                       *
'*          = 7: Formula of Prince-Dormand of 4th/5th order             *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* x        final x-value for iteration. If fehler = 0  we usually have *
'*            x = xend.                                                 *
'* y        final y-values for the solution at x                        *
'* h        final step size used; leave for subsequent calls            *
'* aufrufe  actual number of calls of dgl()                             *
'*                                                                      *
'* Return value (fehler):                                               *
'* =====================                                                *
'* = 0: all ok                                                          *
'* = 1: both error bounds chosen too small for the given mach. constant *
'* = 2: xend <= x0                                                      *
'* = 3: h <= 0                                                          *
'* = 4: n <= 0                                                          *
'* = 5: more right hand side calls than allowed: aufrufe > fmax,        *
'*      x and h contain the current values when stop occured.           *
'* = 6: improper input for embedding formula                            *
'* = 7: lack of available memory (not used here)                        *
'* = 8: Computations completed, but the Prince Dormand formula stiff-   *
'*      ness test indicates possible stiffness.                         *
'* = 9: Computations completed, but both Prince Dormand formula stiff-  *
'*      ness tests indicate possible stiffness. Use method for stiff    *
'*      systems instead '                                               *
'* =10: aufrufe > fmax, see error code 5; AND the Prince Dormand formula*
'*      indicates stiffness; retry using a stiff DE solver '            *
'*                                                                      *
'************************************************************************
Dim MachEps As Double
Dim MACH2 As Double
Dim MACH1 As Double  'machine constant dependent variable which
                      'avoids using too little steps near xend.
Dim xendh As Double  '|xend| - MACH2, carrying same sign as xend
Dim ymax As Double    'Maximum norm of newest approximation of max
                      'order.
Dim hhilf As Double   'aux storage for the latest value of h
                      'produced by step size control. It is saved
                      'here in order to avoid to return a `h' that
                      'resulted from an arbitrary reduction at the
                      'end of the interval.
Dim diff As Double    'distance of the two approximations from the
                      'embedding formula.
Dim s As Double       'indicates acceptance level for results from
                      'embeding formula.

'approximate solution of low order
Dim ybad(n)
'ditto of high order
Dim ygood(n)
   

Dim amEnde As Integer 'flag that shows if the end of the interval
                      'can be reached with the actual step size.
Dim fertig As Integer 'flag indicating end of iterations.

Dim steif1 As Integer 'Flag, that is set in prdo45() if its
                      'stiffness test (dominant eigenvalue)
                      'indicates so. Otherwise no changes.
Dim steifanz As Integer 'counter for number of successive successes
                        'of stiffness test of Shampine and Hiebert in
                        'prdo45().
Dim steif2 As Integer 'Flag, set in prdo45(), when the stiffness
                      'test of  Shampine and Hiebert wa successful
                      'three times in a row; otherwise no changes.
   
'initialize some variables
  fehler = 0
  MachEps = 1.2E-16
  MACH2 = 100 * MachEps
  MACH1 = MachEps ^ 0.75
  amEnde = 0
  fertig = 0
  steif1 = 0
  steif2 = 0
  steifanz = 0
  aufrufe = 1
  ymax = XNormMax(y(), n)

  If xend >= 0# Then
    xendh = xend * (1# - MACH2)
  Else
    xendh = xend * (1# + MACH2)
  End If

' ----------------------- check inputs ----------------------
  If epsabs <= MACH2 * ymax And epsrel <= MACH2 Then
    fehler = 1
    Exit Sub
  End If
  If xendh < x Then
    fehler = 2
    Exit Sub
  End If
  If h < MACH2 * Abs(x) Then
    fehler = 3
    Exit Sub
  End If
  If n <= 0 Then
    fehler = 4
    Exit Sub
  End If
  If methode <> 3 And methode <> 6 And methode <> 7 Then
    fehler = 6
    Exit Sub
  End If
 
' **********************************************************************
' *                                                                    *
' *                       I t e r a t i o n s                          *
' *                                                                    *
' **********************************************************************
  If x + h > xendh Then
                                    'almost at end point ?
    hhilf = h                       'A shortened step might be
    h = xend - x                    'enough.
    amEnde = 1
  End If

  While fertig = 0                  'solve DE system by integrating from
                                    'x0 to xend by suitable steps.
    'choose method
    If methode = 3 Then
      Call ruku23(x, y(), bspn, n, h, ybad(), ygood())
    ElseIf methode = 6 Then
      Call engl45(x, y(), bspn, n, h, ybad(), ygood())
    ElseIf methode = 7 Then
      Call prdo45(x, y(), bspn, n, h, ybad(), ygood(), steif1, steifanz, steif2)
    End If

    aufrufe = aufrufe + methode

    diff = DistMax(ybad(), ygood(), n)

    If (diff < MACH2) Then         'compute s
      s = 2#
    Else
      ymax = XNormMax(ygood(), n)
      s = Sqr(h * (epsabs + epsrel * ymax) / diff)
      If methode <> 3 Then s = Sqr(s)
    End If

    If s > 1# Then                  'integration acceptable ?
      For i = 0 To n - 1            'accept highest order solution
        y(i) = ygood(i)            'move x
      Next i
      
      x = x + h

      If amEnde <> 0 Then           'at end of interval ?
        fertig = 1                  'stop iteration
        If methode = 7 Then
          If steif1 > 0 Or steif2 > 0 Then fehler = 8
          If steif1 > 0 And steif2 > 0 Then fehler = 9
        End If
      ElseIf aufrufe > fmax Then    'too many calls of   dgl() ?
        hhilf = h                   'save actual step size
        fehler = 5                  'report error and stop
        fertig = 1
        If methode = 7 And (steif1 > 0 Or steif2 > 0) Then fehler = 10
      Else                          'Integration was successful
                                    'not at the interval end ?
        h = h * XMin(2#, 0.98 * s)  'increase step size for next
                                    'step properly, at most by
                                    'factor two. Value `0.98*s' is
                                    'used in order to avoid that
                                    'the theoretical value s is
                                    'exceeded by accidental
                                    'rounding errors.
        If x + h > xendh Then      'nearly reached xend ?
          hhilf = h                 '=> One further step with
          h = xend - x              'reduced step size might be
          amEnde = 1                'enough.
          If h < MACH1 * Abs(xend) Then  'very close to xend ?
            fertig = 1                    'finish iteration.
          End If
        End If
      End If
    Else                            'step unsuccessful ?
                                    'before repeating this step
      h = h * XMax(0.5, 0.98 * s)   'reduce step size properly, at
                                    'most by factor 1/2 (for factor
      amEnde = 0                    '0.98: see above).
    End If

  Wend
 
  h = hhilf           'return the latest step size computed by step
                      'size control and  error code to the caller.

End Sub

Function det(n, xmat())
' Determinant  .......................................................
'      integer n                     Dimension of the matrix .........
'      double xmat(0:n-1,0:n-1)       matrix ..........................
 '*====================================================================*
 '*                                                                    *
 '*  det computes the determinant of an n x n real matrix xmat          *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Input parameter:                                                 *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of xmat                                    *
 '*      xmat     matrix(0:n-1,0:n-1)                                  *
 '*               stored in a vector(n*n).                             *
 '*                                                                    *
 '*   Return value:                                                    *
 '*   =============                                                    *
 '*      REAL     Determinant of mat.                                  *
 '*               If the return value = 0, then the matrix is singular *
 '*               or the storage is insufficient                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   subroutine used:                                                 *
 '*   ================                                                 *
 '*                                                                    *
 '*          gaudec ():    LU decomposition of mat                     *
 '*                                                                    *
 '*====================================================================*
 Dim rc As Integer, signd As Integer
 Dim iperm(n), xlu(n * n)
 Dim MachEps As Double
 Dim MAXROOT As Double

  If n < 1 Then   'n not valid
    det = 0#
    Exit Function
  End If
 
  MachEps = 1.2E-16
  MAXROOT = 1E+16

  Call gaudec(n, xmat(), xlu(), iperm(), signd, rc)          'decompose

  If rc <> 0 Or signd = 0 Then
    det = 0#
    Exit Function
  End If

  tmpdet = 1# * signd

  For i = 0 To n - 1
    If Abs(tmpdet) < MachEps Then
      det = 0#
      Exit Function
    ElseIf Abs(tmpdet) > MAXROOT Or Abs(xlu(n * i + i)) > MAXROOT Then
      det = MAXROOT
      Exit Function
    Else
      tmpdet = tmpdet * xlu(n * i + i)                    'compute det
    End If
  Next i
  
  det = tmpdet

End Function

'************************************************************************
'*                                                                      *
'* Test examples for methods to solve first order ordinary systems of   *
'* differential equations.                                              *
'*                                                                      *
'* We supply the right hand sides, their alpha-numeric description and  *
'* the exact analytic solution, if known.                               *
'*                                                                      *
'* When running the main test program, the user can select among the    *
'* examples below.                                                      *
'*                                                                      *
'*                           Quick Basic Release By J-P Moreau, Paris.  *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996".                      *
'************************************************************************
' Note: here, only examples #0 and #3 are implemented in Quick Basic
'
Sub dgl(bspn As Integer, n, x, y(), f())
  If bspn = 0 Then
      f(0) = y(0) * y(1) + Cos(x) - 0.5 * Sin(2# * x)
      f(1) = y(0) ^ 2 + y(1) ^ 2 - (1# + Sin(x))
  ElseIf bspn = 3 Then
      f(0) = y(1)
      f(1) = y(2)
      f(2) = y(3)
      f(3) = y(4)
      f(4) = (45# * y(2) * y(3) * y(4) - 40# * y(3) ^ 3) / (9# * y(2) ^ 2)
  End If
End Sub

'Maximum norm of a difference vector ........
Function DistMax(vector1(), vector2(), n)
'************************************************************************
'* Compute the maximum norm of the difference of two [0..n-1] vectors   *
'*                                                                      *
'* global Name1 used:                                                    *
'* ================                                                     *
'*   None                                                               *
'************************************************************************
' double  abstand     reference value for computation of distance
' double  hilf        distance of two vector elements
  abstand = 0#
  For i = n - 1 To 0 Step -1
    hilf = Abs(vector1(i) - vector2(i))
    If hilf > abstand Then abstand = hilf
  Next i
  DistMax = abstand
End Function

' England's Einbettungs formulas of 4th and 5th degree ................
Sub engl45(x, y(), bspn As Integer, n, h, y4(), y5())
' double    x                   starting point of integration ........
' double    y(0:n-1)            initial value at x ...................
' integer   bspn,               # example
'           n                   number of differential equations .....
' double    h                   step size ............................
' double   y4(0:n-1),           4th order approximation for y at x + h
'          y5(0:n-1)            5th order approximation for y at x + h
' auxiliary vectors
Dim yhilf(n)
Dim k1(n) As Double
Dim k2(n) As Double
Dim k3(n) As Double
Dim k4(n) As Double
Dim k5(n) As Double
Dim k6(n) As Double
'************************************************************************
'* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
'* a solution y at x by using the England embedding formulas on the     *
'* first order system of n differential equations   y' = f(x,y) , as    *
'* supplied by  dgl().                                                  *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    initial x-value                                                 *
'* y    y-values at x, type pVEC                                        *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf, k1..K6: auxiliary vectors                                     *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y4   4th order approximation for y at x + h (pVEC)                   *
'* y5   5th order approximation for y at x + h (pVEC)                   *
'*                                                                      *
'************************************************************************
  Call dgl(bspn, n, x, y(), k1())
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.5 * h * k1(i)
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf(), k2())
  For i = 0 To n - 1
    yhilf(i) = y(i) + (0.25 * h * (k1(i) + k2(i)))
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf(), k3())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h * (-k2(i) + 2# * k3(i))
  Next i
  Call dgl(bspn, n, x + h, yhilf(), k4())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 27# * (7# * k1(i) + 10# * k2(i) + k4(i))
  Next i
  Call dgl(bspn, n, x + (2# / 3#) * h, yhilf(), k5())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 625# * (28# * k1(i) - 125# * k2(i) + 546# * k3(i) + 54# * k4(i) - 378# * k5(i))
  Next i
  Call dgl(bspn, n, x + h / 5#, yhilf(), k6())
  For i = 0 To n - 1
    y4(i) = y(i) + h / 6# * (k1(i) + 4# * k3(i) + k4(i))
    y5(i) = y(i) + h / 336# * (14# * k1(i) + 35# * k4(i) + 162# * k5(i) + 125# * k6(i))
  Next i
End Sub

' Gauss decomposition ................................................
Sub gaudec(n, xmat(), xlumat(), iperm(), signd As Integer, rc As Integer)
' integer   n                        size of matrix ..................
' double    xmat(0:n-1,0:n-1),       Input matrix ....................
'           xlumat(0:n-1,0:n-1)      matrix decomposition ............
' integer   iperm(0:n-1),            row interchanges ................
'           signd,                   sign of perm ....................
'           rc                       error code ......................
'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by
'      vector element [n*i+j].
 '*====================================================================*
 '*                                                                    *
 '*  gaudec decomposes a nonsingular n x n matrix into a product of a  *
 '*  unit lower and an upper triangular matrix. Both triangular factors*
 '*  are stored in lumat (minus the unit diagonal, of course).         *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Input parameter:                                                 *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of  mat and lumat,                         *
 '*               size of  b , x and perm.                             *
 '*      xmat     pointer to original system matrix in vector form.    *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   ==================                                               *
 '*      xlumat   pointer to LU factorization                          *
 '*      iperm    pointer to row permutation vector for xlumat         *
 '*      signd    sign of perm. The determinant of xmat can be computed*
 '*               as the product of the diagonal entries of xlumat     *
 '*               times signd.                                         *
 '*                                                                    *
 '*   Return value (rc):                                               *
 '*   =================                                                *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or invalid input                               *
 '*      = 2      lack of memory                                       *
 '*      = 3      Matrix is singular                                   *
 '*      = 4      Matrix numerically singular                          *
 '*                                                                    *
 '*====================================================================*
  Dim MachEps As Double
  Dim d(n) As Double                  'scaling vector for pivoting

  MachEps = 1.2E-16

  If n < 1 Then
    rc = 1                              'Invalid parameters
    Exit Sub
  End If

  'copy xmat to xlumat
  For i = 0 To n - 1
    For j = 0 To n - 1
      xlumat(n * i + j) = xmat(n * i + j)
    Next j
  Next i

  For i = 0 To n - 1
    iperm(i) = i                         'Initialize iperm
    zmax = 0#
    For j = 0 To n - 1                   'find row maxima
      tmp = Abs(xlumat(n * i + j))
      If tmp > zmax Then zmax = tmp
    Next j

    If zmax = 0# Then                    'xmat is singular
      rc = 3
      Exit Sub
    End If
    d(i) = 1# / zmax
  Next i

  signd = 1                              'initialize sign of iperm

  For i = 0 To n - 1
    piv = Abs(xlumat(n * i + i)) * d(i)
    j0 = i                               'Search for pivot element
    For j = i + 1 To n - 1
      tmp = Abs(xlumat(n * j + i)) * d(j)
      If piv < tmp Then
        piv = tmp                        'Mark pivot element and
        j0 = j                           'its location
      End If
    Next j

    If piv < MachEps Then               'If piv is small, xmat is
      signd = 0                          'nearly singular
      rc = 4
      Exit Sub
    End If

    If j0 <> i Then
      signd = -signd                     'update signd
      Call ISWAP(iperm(j0), iperm(i))    'SWAP1 pivot entries
      Call SWAP1(d(j0), d(i))             'SWAP1 scaling vector
      For j = 0 To n - 1
        'SWAP1 j0-th and i-th rows of xlumat
        Call SWAP1(xlumat(n * j0 + j), xlumat(n * i + j))
      Next j
    End If

    For j = i + 1 To n - 1               'Gauss elimination
      If xlumat(n * j + i) <> 0# Then
        xlumat(n * j + i) = xlumat(n * j + i) / xlumat(n * i + i)
        tmp = xlumat(n * j + i)
        For m = i + 1 To n - 1
          xlumat(n * j + m) = xlumat(n * j + m) - tmp * xlumat(n * i + m)
        Next m
      End If
    Next j
  Next i

  rc = 0  'all ok

End Sub

' Gauss solution .....................................................
Sub gausol(n, xlumat(), iperm(), b(), x(), rc As Integer)
'   integer n                        size of matrix ..................
'   real*8  xlumat(0:n-1,0:n-1)      decomposed matrix (LU) ..........
'   integer perm(0:n-1)              row permutation vector ..........
'   real*8  b(0:n-1),                Right hand side .................
'           x(0:n-1)                 solution ........................
'   integer rc                       error code ......................
 '====================================================================*
 '                                                                    *
 '  gausol  finds the solution x of the linear system  lumat * x = b  *
 '  for the product matrix lumat, that describes an LU decomposition, *
 '  as produced by gaudec.                                            *
 '                                                                    *
 '====================================================================*
 '                                                                    *
 '   Input parameters:                                                *
 '   ================                                                 *
 '      n        integer (n > 0)                                      *
 '               Dimension of xlumat.                                 *
 '      xlumat   Matrix(0:n-1,0:n-1) stored in a vector(n*n).         *
 '               LU factorization, as produced from gaudec.           *
 '      iperm    integer vector(0:n-1)                                *
 '               row permutation vector for xlumat                    *
 '      b        vector(0:n-1)                                        *
 '               Right hand side of the system.                       *
 '                                                                    *
 '   Output parameter:                                                *
 '   ================                                                 *
 '      x        vector(0:n-1)                                        *
 '               Solution vector                                      *
 '                                                                    *
 '   Return value (rc):                                               *
 '   =================                                                *
 '      = 0      all ok                                               *
 '      = 1      n < 1 or other invalid input parameter               *
 '      = 3      improper LU decomposition ( zero diagonal entry)     *
 '                                                                    *
 '====================================================================*
  Dim MachEps As Double
  
  MachEps = 1.2E-16

  If n < 1 Then
     rc = 1                                   'Invalid input parameter
     Exit Sub
  End If

  MachEps = 1.2E-16

  For k = 0 To n - 1                                         'update b
    x(k) = b(iperm(k))
    For j = 0 To k - 1
      x(k) = x(k) - xlumat(n * k + j) * x(j)
    Next j
  Next k

  For k = n - 1 To 0 Step -1                          'back substitute
    Sum = 0#
    For j = k + 1 To n - 1
      Sum = Sum + xlumat(n * k + j) * x(j)
    Next j

    If Abs(xlumat(n * k + k)) < MachEps Then
      rc = 3
      Exit Sub
    End If
    x(k) = (x(k) - Sum) / xlumat(n * k + k)
  Next k

  rc = 0  'all ok

End Sub

' Gauss algorithm for solving linear equations .......................
Sub gauss(mode, n, xmat(), xlumat(), iperm(), b(), x(), signd As Integer, code As Integer)
' integer  mode                      Modus: 0, 1, 2, 3 ...............
' integer  n                         Dimension of matrix .............
' double   xmat(0:n-1,0:n-1),        Input matrix ....................
'          xlumat(0:n-1,0:n-1)       LU decomposition ................
' integer  iperm(0:n-1,0:n-1)        row remutation vector ...........
' double   b(0:n-1),                 right hand side .................
'          x(0:n-1)                  solution of the system ..........
' integer  signd,                    sign of the permutation .........
'          code                      return error code ...............
'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by
'      vector element [n*i+j].
 '*====================================================================*
 '*                                                                    *
 '*  The procedure gauss solves a linear system :  mat * x = b.        *
 '*  Here mat is the nonsingular system matrix, b the right hand side  *
 '*  of the system and x the solution vector.                          *
 '*                                                                    *
 '*  gauss uses the Gauss algorithm and computes a triangular factori- *
 '*  zation of mat and scaled column pivot search.  (Crout method with *
 '*  row SWAP1s).                                                       *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Application:                                                     *
 '*   ============                                                     *
 '*      Solve general linear system with a nonsingular coefficient    *
 '*      matrix.                                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Control parameter:                                               *
 '*   ==================                                               *
 '*      mode     integer;                                             *
 '*               calling modus for gauss:                             *
 '*       = 0     Find factorization and solve linear system           *
 '*       = 1     Find factorization only.                             *
 '*       = 2     Solve linear system only; the factorization is       *
 '*               already available in lumat. This saves work when     *
 '*               solving a linear system repeatedly for several right *
 '*               hand sides and the same system matrix such as when   *
 '*               inverting the matrix.                                *
 '*       = 3     as under 2, additionally we improve the solution     *
 '*               via iterative refinement (not available here).       *
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of mat and lumat,                          *
 '*               size of the vector b, the right hand side, the       *
 '*               solution x and the permutation vector perm.          *
 '*      xmat     matrix of the linear system. It is stored in vector  *
 '*               form.                                                *
 '*      xlumat   (for mode = 2, 3)                                    *
 '*               LU factors of mat                                    *
 '*               xlumat can be stored in the space of xmat.           *
 '*      iperm    (for mode = 2, 3)                                    *
 '*               Permutation vector, of the row interchangfes in      *
 '*               xlumat.                                              *
 '*      b        Right hand side of the system.                       *
 '*      signd    (for mode = 2, 3)                                    *
 '*               sign of the permutation in perm; the determinant of  *
 '*               mat can be computed as the product of the diagonal   *
 '*               entries of lumat times signd.                        *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   ==================                                               *
 '*      xlumat   (for mode = 0, 1)                                    *
 '*               LU factorization of xmat.                            *
 '*      iperm    (for mode = 0, 1)                                    *
 '*               row ermutation vector                                *
 '*      x        (for mode = 0, 2, 3)                                 *
 '*               solution vector(0:n-1).                              *
 '*      signd    (for mode = 0, 1)                                    *
 '*               sign of perm.                                        *
 '*                                                                    *
 '*   Return value (code):                                             *
 '*   ===================                                              *
 '*      =-1      Max. number (MAXITER) of iterative refinements       *
 '*               reached (MAXITER) while mode = 3                     *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or other invalid input                         *
 '*      = 2      lack of memory                                       *
 '*      = 3      Matrix singular                                      *
 '*      = 4      Matrix numerically singular                          *
 '*      = 5      incorrect call                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   subroutines used:                                                *
 '*   ================                                                 *
 '*                                                                    *
 '*      gaudec: determines  LU decomposition                          *
 '*      gausol: solves the linear system                              *
 '*                                                                    *
 '*====================================================================*
Dim rc As Integer

  If (n < 1) Then
    code = 1
    Return
  End If

  'Select mode

    If mode = 0 Then 'Find factorization and solve system ...................
       Call gaudec(n, xmat(), xlumat(), iperm(), signd, rc)
       If rc = 0 Then
         Call gausol(n, xlumat(), iperm(), b(), x(), rc)
         code = rc
       Else
         code = rc
         Exit Sub
       End If

    ElseIf mode = 1 Then 'Find factorization only ...........................
       Call gaudec(n, xmat(), xlumat(), iperm(), signd, rc)
       code = rc
       Exit Sub

    ElseIf mode = 2 Then 'Solve only ........................................
       Call gausol(n, xlumat(), iperm(), b(), x(), rc)
       code = rc
       Exit Sub

    ElseIf mode = 3 Then 'Solve and then use iterative refinement ...........
       Print " fgauss: gausoli not implemented."
       code = 5
       Exit Sub
    End If
  
  code = 5                                               'Wrong call
End Sub

'************************************************************************
'*                                                                      *
'* Solve a first order system of DEs using the implicit Gear method     *
'* of order 4.                                                          *
'*                                                                      *
'* Programming language: ANSI C                                         *
'* Author:               Klaus Niederdrenk, 1.22.1996 (FORTRAN 77)      *
'* Adaptation:           Juergen Dietel, Computing Center, RWTH Aachen  *
'* Source:               FORTRAN 77 source code                         *
'* Date:                 2.26.1996                                      *
'*                                                                      *
'*                       Quick Basic Release By J-P Moreau, Paris.      *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996".                      *
'************************************************************************
'Gear  method of 4th order for DESs of 1st order
Sub gear4(x, xend, bspn As Integer, n, y(), epsabs, epsrel, h, fmax As Integer, aufrufe As Integer, fehler As Integer)
' double x,                      starting or end point ................
' xend                           desired end point (> x) ..............
' integer bspn,                  # of example .........................
' n                              number of DEs ........................
' double y(0:n-1),               initial value or solution ............
' epsabs,                        absolute error bound .................
' epsrel,                        relative error bound .................
' h                              starting or final step size ..........
' integer fmax,                  maximal number of calls of dgl() .....
' aufrufe,                       actual number of calls of dgl() ......
' fehler                         error code ...........................
'************************************************************************
'* Compute the value of the solution at xend, starting with the IVP.    *
'* We use the implicit method of Gear of 4th order with step size       *
'* control which is especially suited for stiff DEs.                    *
'* The local step size control insures that the two error bounds are met*
'* The number of function calls of the right hand side is limited by    *
'* fmax. This function can be used inside a loop to find solutions at   *
'* a specified point to arbitrary accuracy.                             *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x        initial value x0                                            *
'* xend     final value for the integration (xend > x0)                 *
'* bspn     # example                                                   *
'* n        number of DEs                                               *
'* dgl      Function to compute the right hand side f(x0,y0) for the    *
'*          system of DEs (removed from parameters).                    *
'* y        [0..n-1] solution vector y0 of the system of DEs at x0      *
'* epsabs   absolute error bound (>= 0); if zero, we only check for the *
'*          relative error.                                             *
'* epsrel   relative error bound (>= 0); if zero, we only check for the *
'*          absolute error.                                             *
'* h        given step size for first step                              *
'* fmax     maximal number of calls of the right hand side of the system*
'*                                                                      *
'* Output parameters:                                                   *
'* =================                                                    *
'* x        final x-value of the integration; normally equal to xend    *
'* h        final step size; keep for further calls                     *
'* y        [0..n-1]-vector, the solution of the system of DEs at x     *
'* aufrufe  counter of calls of dgl()                                   *
'*                                                                      *
'* erroe code (fehler):                                                 *
'* ====================                                                 *
'*   =   0:  all ok                                                     *
'*   =   1:  Both error bounds too small                                *
'*   =   2:  xend <= x0                                                 *
'*   =   3:  Step size h <= 0                                           *
'*   =   4:  n <= 0                                                     *
'*   =   5:  # calls > fmax; we have not reached the desired xend;      *
'*           repeat function call with actual values of x, y, h.        *
'*   =   6:  Jacobi matrix is singular; x, y, h are the last values     *
'*           that could be computed with accuracy                       *
'*   =   7:  ran out of memory (not used here)                          *
'*   =   8:  error when calling gauss() for the second time;            *
'*           should not occur.                                          *
'*                                                                      *
'************************************************************************
' double eps1,            'MachEps^0.75; used instead of MachEps
'                         'to check whether we have reached xend in
'                         'order to avoid minuscule step sizes
' eps2,                   '100 * MachEps; for comparison with zero
' hs                      'optimal step size for Jacobi matrix
'                         'approximation

Dim hilf(n)               '(0..n-1)-vector
Dim zj(5 * n)             '(0:4,0:n-1)-matrix stored in a vector(5*n)
Dim zjp1(5 * n)           '(0:4,0:n-1)-matrix stored in a vector(5*n)
Dim f(n), ykp1(n)         '(0..n-1)-vectors
Dim fs(n * n), fsg(n * n) '(0..n-1,0..n-1)-matrices stored in vectors(n*n)
Dim con(n)                '(0..n-1)-vector
Dim iperm(n)              '(0..n-1) permutation vector for Gauss elimination

' sg,                      sign of xend
' xe                       |xend| - eps2, carrying the sign of xend
Dim amEnde As Integer     'Flag, indicating that we shall reach
                          'xend with the current step
' ymax,                    Maximum norm of the current
                          'approximate value of y
' dummy,                   aux storage for h
' auxiliary variables (double);
' xka, xke , hka, hk1, diff, eps, q, halt, quot1, quot2, quot3, quot4

Dim done As Integer         'Boolean (0 or 1)
'   nochmal                 'Boolean (0 or 1)
Dim aufrufeawp As Integer
Dim signdet As Integer      'sign of determinant in Gauss

Dim dum(n), dum1(n)       'dummy vectors for gauss
Dim MachEps As Double      'machine smallest real number

'auxiliary vectors
'Dim yhilf(n), k1(n), k2(n), k3(n), k4(n), k5(n), k6(n)

' ------------------------- Initialize ------------------------------
  dummy = 0#: done = 0
  MachEps = 1.2E-16                    'for IBM PC
  eps1 = MachEps ^ 0.75
  eps2 = 100# * MachEps
  hs = 10# * Sqr(MachEps)
  ONE = 1#

  If (xend >= 0#) Then
    sg = 1#
  Else
    sg = -1#
  End If
  xe = (1# - sg * eps2) * xend
  fehler = 0
  aufrufe = 1
  amEnde = 0

' --------- check input parameters -------------------
  ymax = XNormMax(y(), n)
  If epsabs <= eps2 * ymax And epsrel <= eps2 Then
    fehler = 1
  ElseIf xe < x Then
    fehler = 2
  ElseIf h < eps2 * Abs(x) Then
    fehler = 3
  ElseIf n <= 0 Then
    fehler = 4
  End If
  If fehler > 0 Then Exit Sub

' ------------ first integration step ---------------
  If x + h > xe Then
    h = xend - x
    dummy = h
    amEnde = 1
  End If
  For i = 0 To n - 1
    hilf(i) = y(i)
  Next i
  xka = x
  xke = xka
  hka = 0.25 * h
  hk1 = hka

  For k = 1 To 4
    xke = xke + hka

    awp xka, xke, bspn, n, hilf(), epsabs, epsrel, hk1, 6, fmax - aufrufe, aufrufeawp, fehler
        
   aufrufe = aufrufe + aufrufeawp

    If fehler <> 0 Then Exit Sub
    For i = 0 To n - 1
      zjp1(n * k + i) = hilf(i)
    Next i
  Next k

  dgl bspn, n, x, y(), f()
  aufrufe = aufrufe + 1

' ---------- Compute first Gear-Nordsieck approximation -------------------
  For i = 0 To n - 1
    zj(i) = y(i)
    zj(i + n) = h * f(i)
    zj(i + 2 * n) = ONE / 24 * (35 * y(i) - 104 * zjp1(i + n) + 114 * zjp1(i + 2 * n) - 56 * zjp1(i + 3 * n) + 11 * zjp1(i + 4 * n))
    zj(i + 3 * n) = ONE / 12 * (-5 * y(i) + 18 * zjp1(i + n) - 24 * zjp1(i + 2 * n) + 14 * zjp1(i + 3 * n) - 3 * zjp1(i + 4 * n))
    zj(i + 4 * n) = ONE / 24 * (y(i) - 4 * zjp1(i + n) + 6 * zjp1(i + 2 * n) - 4 * zjp1(i + 3 * n) + zjp1(i + 4 * n))
  Next i

' ------------------------ adjust step size --------------------------
  While done = 0
  
    ' --- use Newton method for an implicit approximation ---

    For i = 0 To n - 1
      ykp1(i) = zj(i) + zj(i + n) + zj(i + 2 * n) + zj(i + 3 * n) + zj(i + 4 * n)
    Next i

    dgl bspn, n, x + h, ykp1(), f()

    For k = 0 To n - 1
      'copy vector ykp1 in hilf
      For i = 0 To n - 1
        hilf(i) = ykp1(i)
      Next i
      hilf(k) = hilf(k) - hs
      dgl bspn, n, x + h, hilf(), dum()
      For i = 0 To n - 1
        fs(k * n + i) = dum(i)
      Next i
      For i = 0 To n - 1
        fs(k * n + i) = -h * 0.48 * (f(i) - fs(k * n + i)) / hs
      Next i
      fs(k * n + k) = fs(k * n + k) + ONE
    Next k

    'update number of calls to dgl
    aufrufe = aufrufe + n + 1

    For i = 0 To n - 1
      con(i) = ykp1(i) - 0.48 * (zj(i + n) + 2 * zj(i + 2 * n) + 3 * zj(i + 3 * n) + 4 * zj(i + 4 * n))
      For k = 0 To n - 1
        fsg(k * n + i) = fs(i * n + k)
      Next k
    Next i

    gauss 1, n, fsg(), fsg(), iperm(), dum(), dum1(), signdet, fehler

    If fehler > 0 Then   'error in gauss ?
      fehler = 6
      Exit Sub
    End If

    For iter = 1 To 3
      For i = 0 To n - 1
        hilf(i) = -ykp1(i)
        For k = 0 To n - 1
          hilf(i) = hilf(i) + fs(k * n + i) * ykp1(k)
        Next k
        hilf(i) = h * 0.48 * f(i) + hilf(i) + con(i)
      Next i

      gauss 2, n, fsg(), fsg(), iperm(), hilf(), ykp1(), signdet, fehler

      If fehler > 0 Then
        fehler = 8
        Exit Sub
      End If

      dgl bspn, n, x + h, ykp1(), f()

    Next iter
    'update number of calls to dgl
    aufrufe = aufrufe + 3

'   ---- Compute corresponding Gear-Nordsieck approximation ----

    For i = 0 To n - 1
      hilf(i) = h * f(i) - zj(i + n) - 2 * zj(i + 2 * n) - 3 * zj(i + 3 * n) - 4 * zj(i + 4 * n)
    Next i

    For i = 0 To n - 1
      zjp1(i) = ykp1(i)
      zjp1(i + n) = h * f(i)
      zjp1(i + 2 * n) = zj(i + 2 * n) + 3# * zj(i + 3 * n) + 6# * zj(i + 4 * n) + 0.7 * hilf(i)
      zjp1(i + 3 * n) = zj(i + 3 * n) + 4# * zj(i + 4 * n) + 0.2 * hilf(i)
      zjp1(i + 4 * n) = zj(i + 4 * n) + 0.02 * hilf(i)
    Next i

'   --- decide whether to accept last step ---

'   copy vector zjp1(4) in hilf and zj(4) in con
    For i = 0 To n - 1
      hilf(i) = zjp1(i + 4 * n)
      con(i) = zj(i + 4 * n)
    Next i

    diff = DistMax(hilf(), con(), n)
    ymax = XNormMax(ykp1(), n)
    eps = (epsabs + epsrel * ymax) / 6#
    q = Sqr(Sqr(eps / diff)) / 1.2

    If (diff < eps) Then

      ' --- accept last step; prepare for next one ---

      x = x + h
      'copy vector ykp1 in y
      For i = 0 To n - 1
        y(i) = ykp1(i)
      Next i

'     stop integration, if interval end xend has been reached
'     or if there has been too many function dgl calls.

      nochmal = 0
      While nochmal = 0
        If amEnde <> 0 Then
          h = dummy
          Exit Sub
        ElseIf aufrufe > fmax Then
          fehler = 5
          Exit Sub
        End If

        ' --- adjust step size for next step ---
        halt = h
        h = XMin(q, 2#) * h

        If x + h >= xe Then
          dummy = h
          h = xend - x
          amEnde = 1

          ' --- close enough to xend => stop integration ---
          If h < eps1 * Abs(xend) Then nochmal = 1
        End If
        If nochmal = 0 Then GoTo 10
      Wend

'     ------ compute Gear-Nordsieck approximation -----
'     ------ for the next step                    -----
10    quot1 = h / halt
      quot2 = quot1 ^ 2
      quot3 = quot2 * quot1
      quot4 = quot3 * quot1
      For i = 0 To n - 1
        zj(i) = zjp1(i)
        zj(i + n) = quot1 * zjp1(i + n)
        zj(i + 2 * n) = quot2 * zjp1(i + 2 * n)
        zj(i + 3 * n) = quot3 * zjp1(i + 3 * n)
        zj(i + 4 * n) = quot4 * zjp1(i + 4 * n)
      Next i
    Else
'     ------ repeat last step with smaller step size;  -----
'     -------- adjust Gear-Nordsieck approximation ---------
      halt = h
      h = XMax(0.5, q) * h
      quot1 = h / halt
      quot2 = quot1 ^ 2
      quot3 = quot2 * quot1
      quot4 = quot3 * quot1
      For i = 0 To n - 1
        zj(i + n) = quot1 * zj(i + n)
        zj(i + 2 * n) = quot2 * zj(i + 2 * n)
        zj(i + 3 * n) = quot3 * zj(i + 3 * n)
        zj(i + 4 * n) = quot4 * zj(i + 4 * n)
      Next i
      amEnde = 0
    End If
  Wend  'while done=0

End Sub 'gear4

Sub ISWAP(ia, ib)
Dim tnp As Integer
  tmp = ib: ib = ia: ia = tmp
End Sub

' Gauss for multiple right hand sides ................................
Sub mgauss(n, k, xmat(), rmat(), code As Integer)
'   integer n,                       Dimension of system .............
'           k                        number of right hand sides ......
'   real*8  mat(0:n-1,0:n-1),        original matrix .................
'           rmat(0:n-1,0:n-1)        Right hand sides/solutions ......
'   integer code                     Error code ......................
 '*====================================================================*
 '*                                                                    *
 '*  mgauss  finds the solution matrix x for the linear system         *
 '*  mat * x = rmat with an  n x n coefficient matrix mat and a        *
 '*  n x k matrix rmat of right hand sides. Here mat must be           *
 '*  nonsingular.                                                      *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of xmat.                                   *
 '*      k        integer k;  (k > 0)                                  *
 '*               number of right hand sides                           *
 '*      xmat     matrix(0:n-1,0:n-1) stored in a vector(n*n)          *
 '*               n x n original system matrix                         *
 '*      rmat     matrix(0:n-1,0:n-1) stored in a vector(n*n)          *
 '*               Right hand sides                                     *
 '*                                                                    *
 '*   Output parameter:                                                *
 '*   ================                                                 *
 '*      rmat     solution matrix for the system.                      *
 '*               The input right hand sides are lost.                 *
 '*                                                                    *
 '*   Return value (code):                                             *
 '*   ===================                                              *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or k < 1 or invalid input parameter            *
 '*      = 2      lack of memory (not used here)                       *
 '*      = 3      mat is numerically singular.                         *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   subroutine used:                                                 *
 '*   ================                                                 *
 '*                                                                    *
 '*      gaudec:  LU decomposition of mat.                             *
 '*                                                                    *
 '*====================================================================*
 Dim signd As Integer, rc As Integer
 Dim iperm(n)
 Dim xlu(n * n), x(n)
 Dim MachEps As Double
 
  If n < 1 Or k < 1 Then                        'Invalid parameter
    code = 1
    Exit Sub
  End If

  MachEps = 1.2E-16

  Call gaudec(n, xmat(), xlu(), iperm(), signd, rc) 'compute factorization
                                                    'in matrix lu
  If rc <> 0 Or signd = 0 Then                      'if not possible
    code = 3                                        'exit with code=3
    Exit Sub
  End If

  For m = 0 To k - 1                   'Loop over the right hand sides
    For i = 0 To n - 1                               'Updating the b's
      x(i) = rmat(n * iperm(i) + m)
      For j = 0 To i - 1
        x(i) = x(i) - xlu(n * i + j) * x(j)
      Next j
    Next i

    For i = n - 1 To 0 Step -1                      'back substitution
      Sum = 0#
      For j = i + 1 To n - 1
        Sum = Sum + xlu(n * i + j) * x(j)
      Next j

      If Abs(xlu(n * i + i)) < MachEps Then  'invalid LU decomposition
        code = 2
        Exit Sub
      End If
      x(i) = (x(i) - Sum) / xlu(n * i + i)
    Next i

    For j = 0 To n - 1                                    'Save result
      rmat(n * j + m) = x(j)
    Next j
  Next m

  code = 0

End Sub

' embedding formulas of Prince-Dormand of 4./5. order .........................
Sub prdo45(x, y(), bspn As Integer, n, h, y4(), y5(), steif1 As Integer, steifanz As Integer, steif2 As Integer)
' double    x,                      starting point of integration .....
'           y(0:n-1)                initial value at x ................
' integer   bspn,                   # example .........................
'           n                       number of DEs .....................
' double    h,                      step size .........................
'           y4(0:n-1),              solution of 4th order at x+h ......
'           y5(0:n-1)               solution of 5th order at x+h ......
' auxiliary flags
' integer   steif1, steifanz,steif2
' auxiliary vectors
Dim yhilf(n)
Dim k1(n) As Double, k2(n) As Double
Dim k3(n) As Double, k4(n) As Double
Dim k5(n) As Double, k6(n) As Double
Dim k7(n) As Double
Dim g6(n) As Double, g7(n) As Double
'************************************************************************
'* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
'* a solution y at x by using the Prince-Dormand embedding formulas on  *
'* the first order system of n differential equations y' = f(x,y) , as  *
'* supplied by  dgl().                                                  *
'* Simultaneously we perform two tests for stiffness whose results are  *
'* stored in steif1 and steif2.                                         *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    initial x-value                                                 *
'* y    y-values at x (pVEC)                                            *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf, k1..k7,g6,g7: auxiliary vectors.                              *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y4   4th order approximation for y at x + h                          *
'* y5   5th order approximation for y at x + h                          *
'*                                                                      *
'***********************************************************************}
  Dim steifa As Integer  'Flag which is set if the second test for stiffness
                         'Shampine und Hiebert) is positive; otherwise the
                         'flag is erased.

  Call dgl(bspn, n, x, y(), k1())
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.2 * h * k1(i)
  Next i
  Call dgl(bspn, n, x + 0.2 * h, yhilf(), k2())
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.075 * h * (k1(i) + 3# * k2(i))
  Next i
  Call dgl(bspn, n, x + 0.3 * h, yhilf(), k3())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 45# * (44# * k1(i) - 168# * k2(i) + 160# * k3(i))
  Next i
  Call dgl(bspn, n, x + 0.8 * h, yhilf(), k4())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 6561# * (19372# * k1(i) - 76080# * k2(i) + 64448# * k3(i) - 1908# * k4(i))
  Next i
  Call dgl(bspn, n, x + (8# / 9#) * h, yhilf(), k5())
  For i = 0 To n - 1
    g6(i) = y(i) + h / 167904# * (477901# * k1(i) - 1806240# * k2(i) + 1495424# * k3(i) + 46746# * k4(i) - 45927# * k5(i))
  Next i
  Call dgl(bspn, n, x + h, g6(), k6())
  For i = 0 To n - 1
    g7(i) = y(i) + h / 142464# * (12985# * k1(i) + 64000# * k3(i) + 92750# * k4(i) - 45927# * k5(i) + 18656# * k6(i))
  Next i
  Call dgl(bspn, n, x + h, g7(), k7())
  For i = 0 To n - 1
    y5(i) = g7(i)
    y4(i) = y(i) + h / 21369600# * (1921409# * k1(i) + 969088# * k3(i) + 13122270# * k4(i) - 5802111# * k5(i) + 1902912# * k6(i) + 534240# * k7(i))
  Next i

' Test for stiffness via dominant eigenvalue

  If DistMax(k7(), k6(), n) > 3.3 * DistMax(g7(), g6(), n) Then steif1 = 1

' one step in steffness test of Shampine & Hiebert

  For i = 0 To n - 1
    g6(i) = h * (2.2 * k2(i) + 0.13 * k4(i) + 0.144 * k5(i))
    g7(i) = h * (2.134 * k1(i) + 0.24 * k3(i) + 0.1 * k6(i))
  Next i

  If DistMax(g6(), g7(), n) < DistMax(y4(), y5(), n) Then
    steifa = 1
  Else
    steifa = 0
  End If

  If (steifa > 0) Then
    steifanz = steifanz + 1
    If steifanz >= 3 Then steif2 = 1
  Else
    steifanz = 0
  End If

End Sub

'print a REAL square matrix with Name1 (debug only)
Sub PrintMat(Name1 As String, n, xmat())
' Name1       Matrix caption
' n           Size of matrix
' double xmat(0:n-1,0:n-1) stored in a vector(n*n)
' matrix to be printed
  Print " "; Name1
  For i = 0 To n - 1
    For j = 0 To n - 1
      Print " "; xmat(n * i + j);
    Next j
    Print
  Next i
End Sub

'debug only
Sub PrintVec(Name1 As String, n, V())
  Print Name1
  For i = 0 To n - 1
    Print " "; V(i);
  Next i
  Print
End Sub

' Runge-Kutta embedding formulas of 2nd, 3rd degree ....................
Sub ruku23(x, y(), bspn As Integer, n, h, y2(), y3())
'************************************************************************
'* Compute 2nd and 3rd order approximates y2, y3 at x + h starting with *
'* a solution y at x by using Runge-Kutta embedding formulas on the     *
'* first order system of n differential equations   y' = f(x,y) , as    *
'* supplied by  dgl().                                                  *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    x-value of left end point                                       *
'* y    y-values at x                                                   *
'* bspn # example                                                       *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf,k1,k2,k3: auxiliary vectors defined in module awp.             *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y2   2nd order approximation for y at x + h                          *
'* y3   3rd order approximation for y at x + h                          *
'*                                                                      *
'************************************************************************
  Dim yhilf(n)
  Dim k1(n) As Double
  Dim k2(n) As Double
  Dim k3(n) As Double
  Call dgl(bspn, n, x, y(), k1())
  For i = 0 To n - 1
    yhilf(i) = y(i) + h * k1(i)
  Next i
  Call dgl(bspn, n, x + h, yhilf(), k2())
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.25 * h * (k1(i) + k2(i))
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf(), k3())
  For i = 0 To n - 1
    y2(i) = y(i) + 0.5 * h * (k1(i) + k2(i))
    y3(i) = y(i) + h / 6# * (k1(i) + k2(i) + 4# * k3(i))
  Next i
End Sub

'***************************************************************
'* This subroutine defines n, number of DEs and dgltxt(), the  *
'* text decription of current example.                         *
'***************************************************************
'NOTE: here, only examples #0, #3 are implemented.
  Sub settxt(bspn As Integer, n, dgltxt$())
    If bspn = 0 Then
      n = 2
      dgltxt$(0) = " y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)"
      dgltxt$(1) = " y2' = y1 * y1 + y2 * y2 - (1 + sin(x))"
    ElseIf bspn = 3 Then
      n = 5
      dgltxt$(0) = " y1' = y2"
      dgltxt$(1) = " y2' = y3"
      dgltxt$(2) = " y3' = y4"
      dgltxt$(3) = " y4' = y5"
      dgltxt$(4) = " y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)"
    End If
  End Sub

Sub SWAP1(a, b)
  tmp = b: b = a: a = tmp
End Sub

Function XMax(a, b)
  If a >= b Then
    XMax = a
  Else
    XMax = b
  End If
End Function

Function XMin(a, b)
  If a <= b Then
    XMin = a
  Else
    XMin = b
  End If
End Function

' Find the maximum norm of a REAL vector ............................
Function XNormMax(vektor(), n)
'      double  vektor(0:n-1)                 vector .................
'      integer n                             length of vector .......
' ************************************************************************
' *          Return the maximum norm of a [0..n-1] vector v              *
' *                                                                      *
' ************************************************************************
Dim norm As Double                            ' local max.
'   double betrag                             ' magnitude of a component
    norm = 0#
    For i = 0 To n - 1
      betrag = Abs(vektor(i))
      If betrag > norm Then norm = betrag
    Next i
    XNormMax = norm
End Function

' -------------------------- END  mgear1.bas -------------------------