DECLARE FUNCTION XNormMax# (vektor#(), n%)
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%)
'*******************************************************************
'* Solve an ordinary system of first order differential equations  *
'* using subroutine AWP with step size and precision controls.     *
'* --------------------------------------------------------------- *
'* SAMPLE RUNS:                                                    *
'* Example #1                                                      *
'*                                                                 *
'* # example (0 to 3) ..........: 0                                *
'* Number of equations in system: 2                                *
'* Starting x ..................: 0                                *
'* Ending x ....................: 2                                *
'* Absolute precision ..........: 1D-8                             *
'* Relative precision ..........: 1D-10                            *
'* Starting integration step ...: 0.1                              *
'* Max. number of function calls: 1000                             *
'*                                                                 *
'* Y(0) = 0                                                        *
'* Y(1) = 0                                                        *
'*                                                                 *
'* D.E. System to solve:                                           *
'* y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)                       *
'* y2' = y1 * y1 + y2 * y2 - (1 + sin(x))                          *
'*                                                                 *
'* x= 2                                                            *
'* Y 1 =  9.193162467174124D-02                                    *
'* Y 2 = -1.363855036258106                                        *
'*                                                                 *
'* Number of function calls:  685                                  *
'* Final step size ........:  1.728309648562938D-02                *
'* Error code .............:  0                                    *
'*                                                                 *
'* Example #2                                                      *
'*                                                                 *
'* # example (0 to 3) ..........: 1                                *
'* Number of equations in system: 1                                *
'* Starting x ..................: 0                                *
'* Ending x ....................: 1                                *
'* Absolute precision ..........: 1D-8                             *
'* Relative precision ..........: 1D-10                            *
'* Starting integration step ...: 0.1                              *
'* Max. number of function calls: 500                              *
'*                                                                 *
'* Y(0) = 1                                                        *
'*                                                                 *
'* D.E. System to solve:                                           *
'* y' = -y + x/((1+x)*(1+x))                                       *
'*                                                                 *
'* x=  1                                                           *
'* Y(0) =  .4999999999975014                                       *
'*                                                                 *
'* Number of function calls:  187                                  *
'* Final step size ........:  4.584567179138789D-02                *
'* Error code .............:  0                                    *
'*                                                                 *
'* Example #3                                                      *
'*                                                                 *
'* # example (0 to 3) ..........: 3                                *
'* Number of equations in system: 5                                *
'* Starting x ..................: 0                                *
'* Ending x ....................: 2                                *
'* Absolute precision ..........: 1D-8                             *
'* Relative precision ..........: 1D-10                            *
'* Starting integration step ...: 0.1                              *
'* Max. number of function calls: 2000                             *
'*                                                                 *
'* Y(0) = 1                                                        *
'* Y(1) = 1                                                        *
'* Y(2) = 1                                                        *
'* Y(3) = 1                                                        *
'* Y(4) = 1                                                        *
'*                                                                 *
'* D.E. System to solve:                                           *
'* y1' = y2                                                        *
'* y2' = y3                                                        *
'* y3' = y4                                                        *
'* y4' = y5                                                        *
'* y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)   *
'*                                                                 *
'* x=  2                                                           *
'* Y(0) =  6.708203932537073                                       *
'* Y(1) =  5.341640786529041                                       *
'* Y(2) =  2.414953415673884                                       *
'* Y(3) = -1.448972049484034                                       *
'* Y(4) = -1.448972049312524                                       *
'*                                                                 *
'* Number of function calls:  1387                                 *
'* Final step size ........:  9.544375718212773D-03                *
'* Error code .............:  0                                    *
'*                                                                 *
'* --------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges   *
'*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].      *
'*                                                                 *
'*                  Quick Basic Release 1.1 By J-P Moreau, Paris.  *
'*                                 (www.jpmoreau.fr)               *
'*                                                                 *
'* Release 1.1 (12/21/06): added two more examples.                *
'*******************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

DIM fmax AS INTEGER
DIM bspnum AS INTEGER


  CLS
  PRINT
  INPUT " # Example ...............: ", bspnum   'see sub dgl
  INPUT " Number of DEs ...........: ", n        'number of equations in system
  DIM Y(n)
  INPUT " Starting x value ........: ", xbegin
  INPUT " Ending x value ..........: ", xend
  INPUT " Absolute precision ......: ", epsabs
  INPUT " Relative precision ......: ", epsrel
  INPUT " Starting integration step: ", h
  INPUT " Max. number of dgl calls : ", fmax
  FOR i = 0 TO n - 1                             'initial value of Y
    PRINT " Y("; i; ") = "; : INPUT "", Y(i)
  NEXT i
  
  method = 6        'England embedding formulas

  'call integration procedure with step size control
  awp xbegin, xend, bspnum, n, Y(), epsabs, epsrel, h, method, fmax, ncalls, ierror

  'write results
  CLS
  PRINT
  PRINT " x= "; xend
  FOR i = 0 TO n - 1
    PRINT " y"; i + 1; " = "; Y(i)
  NEXT i
  PRINT
  PRINT " Number of function calls: "; ncalls
  PRINT " Final step size ........: "; h
  PRINT " Error code .............: "; ierror
  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 ^ .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#, .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(.5, .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

'************************************************************************
'*                                                                      *
'* 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, #1 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) - .5 * SIN(2# * x)
      f(1) = Y(0) ^ 2 + Y(1) ^ 2 - (1# + SIN(x))
  ELSEIF bspn = 1 THEN
      f(0) = -Y(0) + x / ((1# + x) * (1# + 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) + .5 * h * k1(i)
  NEXT i
  CALL dgl(bspn, n, x + .5 * h, yhilf(), k2())
  FOR i = 0 TO n - 1
    yhilf(i) = Y(i) + (.25 * h * (k1(i) + k2(i)))
  NEXT i
  CALL dgl(bspn, n, x + .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

SUB ISWAP (ia, ib)
DIM tnp AS INTEGER
  tmp = ib: ib = ia: ia = tmp
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) + .2 * h * k1(i)
  NEXT i
  CALL dgl(bspn, n, x + .2 * h, yhilf(), k2())
  FOR i = 0 TO n - 1
    yhilf(i) = Y(i) + .075 * h * (k1(i) + 3# * k2(i))
  NEXT i
  CALL dgl(bspn, n, x + .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 + .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) + .13 * k4(i) + .144 * k5(i))
    g7(i) = h * (2.134 * k1(i) + .24 * k3(i) + .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) + .25 * h * (k1(i) + k2(i))
  NEXT i
  CALL dgl(bspn, n, x + .5 * h, yhilf(), k3())
  FOR i = 0 TO n - 1
    y2(i) = Y(i) + .5 * h * (k1(i) + k2(i))
    y3(i) = Y(i) + h / 6# * (k1(i) + k2(i) + 4# * k3(i))
  NEXT i
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