DECLARE SUB dgl (bspn AS INTEGER, n%, x#, y#(), f#())
DECLARE SUB extrapol (row AS INTEGER, fhilf#(), y#(), n%, epsabs#, epsrel#, x#, h#, h0#, ahead AS INTEGER, index%, extmax AS INTEGER, bufol() AS INTEGER)
DECLARE FUNCTION ibulstoe% (num%, x#, xend#, n%, y#(), epsabs#, epsrel#, h#, hmax#, neu%, fmax AS INTEGER, aufrufe AS INTEGER)
DECLARE FUNCTION Min% (ia%, ib%)
DECLARE FUNCTION XMax# (A#, b#)
DECLARE FUNCTION XMin# (A#, b#)
'*******************************************************************
'* Solve an ordinary system of first order differential equations  *
'* using the extrapolation method of Bulirsch-Stoer-Gragg.         *
'* --------------------------------------------------------------- *
'* 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.05                             *
'* 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 =  .0919316246481357                                        *
'* Y 2 = -1.36385503607559                                         *
'*                                                                 *
'* Number of function calls:  578                                  *
'* Final step size ........:  2.522919442059587D-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 1 =  .5000000001223267                                        *
'*                                                                 *
'* Number of function calls:  192                                  *
'* Final step size ........:  2.819970295498907D-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 1 =  6.708203932363743                                        *
'* Y 2 =  5.341640786173861                                        *
'* Y 3 =  2.414953415257717                                        *
'* Y 4 = -1.448972049292256                                        *
'* Y 5 = -1.448972048644391                                        *
'*                                                                 *
'* Number of function calls:  930                                  *
'* Final step size ........:  3.768565498571516D-02                *
'* 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.                *
'*******************************************************************
'Program Test_Bulirsch
DEFINT I-N
DEFDBL A-H, O-Z

OPTION BASE 0

  DIM erreur AS INTEGER
  DIM bspnum AS INTEGER
  DIM fmax 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

  'call integration procedure with the extrapolation method
  erreur = ibulstoe(bspnum, xbegin, xend, n, y(), epsabs, epsrel, h, 1#, 1, fmax, ncalls)

  '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 .............: "; erreur

END 'of main program

'Utility subroutine
SUB copyvector (A(), b(), n)
  FOR i = 0 TO n - 1
    A(i) = b(i)
  NEXT i
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

SUB extrapol (row AS INTEGER, fhilf(), y(), n, epsabs, epsrel, x, h, h0, ahead AS INTEGER, index, extmax AS INTEGER, bufol() AS INTEGER)
'    row: integer
'    fhilf(12,2)
'    y(n)
'    ahead: boolean
'    bufol(12): integer (bulirsch sequence)
'************************************************************************
'* Perform one extrapolation step for function ibulstoe                 *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* row     pointer to the arrays  bufol and fhilf                       *
'* fhilf   pointer to Matrix for extrapolation values (here vector)     *
'* y       y-value of solution at x                                     *
'* n       number of equations in the system                            *
'* epsabs  absolute error bound                                         *
'* epsrel  relative error bound                                         *
'* x       x-value                                                      *
'* h       local step sizes                                             *
'* h0                                                                   *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* ahead  *ahead  = 1:  Step is acceptable                              *
'*        *ahead != 1:  Step not acceptable                             *
'* index  pointer to the arrays bufol and fhilf                         *
'* x      final x-value                                                 *
'* h      final step size                                               *
'*                                                                      *
'* global names used:                                                   *
'* ==================                                                   *
'*        extmax, bufol, XMin, XMax, Copyvector.                        *
'************************************************************************
DIM column AS INTEGER  'column of extrapolation scheme
'   diffmax,            maximal difference of two columns in scheme
'   fhilfi1,            column in extrapolation scheme
'   fhilfi2,            column to the right of fhilfi1
'   bufolr,             element in Bulirsch sequence
'   bufoli,             idem
'   ymax,               maximal value in a column
'   help                auxiliary variable

  ymax = 0#
  maxcol = Min(row + 1, extmax)
  ahead = 0  'False

  FOR column = 2 TO maxcol
    IF ahead <> 0 THEN EXIT SUB
    index = Min(11 - column, row - column + 3)
    diffmax = 0#
    FOR i = 0 TO n - 1
      fhilfi1 = fhilf((index - 1) * n + i)
      fhilfi2 = fhilf((index - 2) * n + i)
      bufolr = 1# * bufol(row)
      bufoli = 1# * bufol(index - 2)
      fhilf((index - 2) * n + i) = fhilfi1 + (fhilfi1 - fhilfi2) / ((bufolr / bufoli) ^ 2 - 1#)
      fhilfi1 = fhilf((index - 1) * n + i)
      fhilfi2 = fhilf((index - 2) * n + i)
      ymax = XMax(ymax, ABS(fhilfi2))
      diffmax = XMax(diffmax, ABS(fhilfi2 - fhilfi1))
    NEXT i
   
    IF diffmax < epsrel * ymax + epsabs THEN    'Step acceptable ?
      x = x + h0
      'Copyvector(y, fhilf(index - 2, n)
      FOR i = 0 TO n - 1
        y(i) = fhilf((index - 2) * n + i)
      NEXT i
      help = 1# * (column - extmax)
      'Step size for next step
      h = .9# * h * (.6 ^ help)
      row = -1
      ahead = 1   'True
    END IF
  NEXT column
END SUB

'Extrapolation method for 1st order DE systems.
FUNCTION ibulstoe (num, x, xend, n, y(), epsabs, epsrel, h, hmax, neu, fmax AS INTEGER, aufrufe AS INTEGER)
'        num: integer            # example in Sub dgl ...............
'        x: Double               initial x-value/ final x-value .....
'        xend: Double            desired end point ..................
'        n: integer              number of DEs ......................
'        y(n)                    initial y-value/ final y-value .....
'        epsabs,                 absolute error bound ...............
'        epsrel: Double          relative error bound ...............
'        h: Double               initial/final step size ............
'        hmax: Double            maximal step size ..................
'        neu: boolean            use an outside given x ? ...........
'        fmax: integer           maximal # of calls of  dgl .........
'        aufrufe: integer        actual # of calls of dgl ...........
'************************************************************************
'* Compute the solution of an intial value problem for a first order    *
'* system of ordinary differential equations                            *
'*                                                                      *
'*      y' = f(x,y)        with initial value y(x0) = y0                *
'*                                                                      *
'* by using the extrapolation method of Bulirsch-Stoer-Gragg.           *
'* The maximal extrapolation order is set internally depending on the   *
'* machine constant; the step size control follows the ideas in :       *
'*     HALL, G.; WATT, J.M: Modern Numerical Methods for Ordinary       *
'*                          Differential Equations, Oxford 1976,        *
'*                          Clarendon Press,  [HALL76]                  *
'* REMARK :                                                             *
'* The dynamic allocations from this function are only set free for a   *
'* new run or with improper input (fehler = TRUE or Return value >= 2). *
'* Hence we advise the user to end work with this function by one call  *
'* with an improper value such as n = -2 in order to free up all storage*
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x        x-value                                                     *
'* y        [0..n-1] vector with initial y-value at x                   *
'* n        number of equations in the system                           *
'* xend     desired final value for x, xend < x is allowed              *
'* h        step size for next step                                     *
'* hmax     maximal step size; hmax > 0                                 *
'* epsabs\  error bounds for absolute and relative errors, both non     *
'* epsrel/  negative.                                                   *
'*          We test for the mixed error:                                *
'*              |lokaler Fehler|  <=  |y| * epsrel + epsabs.            *
'*          If epsrel = 0, we test the absolute error; if epsabs = 0 we *
'*          test the relative error. epsrel will be set to 10 * machine *
'*          constant if it has been chosen smaller than this.           *
'* neu      Flag, indicating whether this function is called to compute *
'*          on after one partial run that was cut short due to too many *
'*          evaluations:                                                *
'*          neu = FALSE: continue with old values                       *
'*          neu = TRUE:  start new                                      *
'* fmax     upper bound for the number of allowed function evaluations  *
'*          of the right hand side via dgl()                            *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* x        final x-value of integration; usually x = xend              *
'* y        [0..n-1] y-value vector at x                                *
'* h        final step size; should be used unchanged in the next call  *
'*          if neu = TRUE; should be redefined for neu = TRUE.          *
'* aufrufe  number of calls of dgl()                                    *
'*                                                                      *
'* Return value :                                                       *
'* ==============                                                       *
'* = 0: all ok; after a new setting for xend, the function can be used  *
'*      again for the same problem, or all inputs may be changed        *
'* = 1: method has not reached xend after  FMAX right hand side         *
'*      evaluations; try another call with same parameters or change    *
'*      error bounds.                                                   *
'* = 2: step size less than 4 * machine constant.                       *
'*      For a future call, increase step size and error bounds.         *
'* = 3: epsabs or epsrel negative, or both zero.                        *
'* = 4: xend = x                                                        *
'* = 5: hmax negative.                                                  *
'* = 6: n <= 0.                                                         *
'*                                                                      *
'* global names used:                                                   *
'* ==================                                                   *
'* extrapol, MACHEPS, XMin, Min, XMax, copyvector, ABS.                 *
'************************************************************************
DIM MACHEPS AS DOUBLE   ' Smallest machine-dependant real

DIM bufol(12)  AS INTEGER 'The Bulirsch sequence
DIM fhilf(12 * n)       ' Index 0..7: Matrix for the extrapolation scheme
                        ' Index 8..11: aux vectors
'    x0 as Double         x-value on call of dgl
DIM row  AS INTEGER     ' row of extrapolation scheme
'    absh,                magnitude of step size
'    absh0,               magnitude of aux step size h0
'    h0,                  auxiliary step size
'    hilf: Double         auxiliary variable
DIM count  AS INTEGER   ' Loop variable for mid-point rule
'    index: integer       Index for the vectors fhilf and bufol
DIM ahead  AS INTEGER   ' Flag indicating acceptance of step
DIM extmax AS INTEGER   ' Maximum order of extrapolation

DIM tmp(n), tmp1(n)     ' Auxiliary vectors of size n for calls of dgl

  MACHEPS = 2.225D-16

  'Initialize bufol array
  bufol(0) = 2: bufol(1) = 4: bufol(2) = 6: bufol(3) = 8: bufol(4) = 12
  bufol(5) = 16: bufol(6) = 24: bufol(7) = 32: bufol(8) = 48
  bufol(9) = 64: bufol(10) = 96: bufol(11) = 128

  FOR i = 0 TO 11
    FOR j = 0 TO n - 1
      fhilf(i * n + j) = 0#
    NEXT j
  NEXT i

  h0 = 0#

  'check input
  IF epsabs < 0# OR epsrel < 0# OR epsabs + epsrel <= 0# THEN
    ibulstoe = 3
    EXIT FUNCTION
  ELSEIF xend = x THEN
    ibulstoe = 4
    EXIT FUNCTION
  ELSEIF hmax <= 0# THEN
    ibulstoe = 5
    EXIT FUNCTION
  END IF
  IF n <= 0 THEN
    ibulstoe = 6
    EXIT FUNCTION
  END IF

  IF neu = 0 THEN 'new call or repeat call due to excessive evaluations ?
    x = fhilf(9 * n)
    'copyvector(y, fhilf(8), n)
    FOR i = 0 TO n - 1
      y(i) = fhilf(8 * n + i)
    NEXT i
  END
  ELSE
    row = -1
  END IF
 
  aufrufe = 0
  ahead = 1    'TRUE
  extmax = INT(-LOG(MACHEPS) / LOG(2#) / 7.5# + .5#)

  epsrel = XMax(epsrel, 10# * MACHEPS)

  'main integration loop
15  IF neu <> 0 THEN                                      'new calculation

      IF ahead <> 0 THEN                                       'new step ?
        absh = ABS(h)
        absh = XMin(absh, hmax)
        h0 = XMin(absh, ABS(xend - x))
        IF xend < x THEN h0 = -h0
        absh0 = ABS(h0)

        IF absh0 <= 4# * MACHEPS * ABS(x) THEN
          ibulstoe = 0
          EXIT FUNCTION                                       'normal exit
        END IF
        ahead = 0  'FALSE
      END IF


20    row = row + 1                       'find step size for extrapolation
        h = h0 / bufol(row)
        x0 = x                            'Euler step; save initial values

        'copyvector(fhilf(8), y, n)
        FOR i = 0 TO n - 1
          fhilf(8 * n + i) = y(i)
        NEXT i

        FOR i = 0 TO n - 1
          tmp(i) = fhilf(8 * n + i)
        NEXT i

        dgl num, n, x0, tmp(), tmp1()

        FOR i = 0 TO n - 1
          fhilf(11 * n + i) = tmp1(i)
        NEXT i

        FOR i = 0 TO n - 1
          fhilf(9 * n + i) = fhilf(8 * n + i) + h * fhilf(11 * n + i)
        NEXT i

        x0 = x0 + h

        FOR i = 0 TO n - 1
          tmp(i) = fhilf(9 * n + i)
        NEXT i

        dgl num, n, x0, tmp(), tmp1()

        FOR i = 0 TO n - 1
          fhilf(11 * n + i) = tmp1(i)
        NEXT i

        'use mid-point rule
        FOR count = 1 TO bufol(row) - 1
        
          FOR i = 0 TO n - 1
            fhilf(10 * n + i) = fhilf(8 * n + i) + 2# * h * fhilf(11 * n + i)
          NEXT i

          x0 = x0 + h

          FOR i = 0 TO n - 1
            tmp(i) = fhilf(10 * n + i)
          NEXT i

          dgl num, n, x0, tmp(), tmp1()

          FOR i = 0 TO n - 1
            fhilf(11 * n + i) = tmp1(i)
          NEXT i

          FOR j = 8 TO 10                         'store for next step
            'copyvector(fhilf(j), fhilf(j + 1), n)
            FOR i = 0 TO n - 1
              fhilf(j * n + i) = fhilf((j + 1) * n + i)
            NEXT i
          NEXT j

        NEXT count

        FOR i = 0 TO n - 1
          tmp(i) = fhilf(9 * n + i)
        NEXT i

        dgl num, n, x0, tmp(), tmp1()

        FOR i = 0 TO n - 1
          fhilf(11 * n + i) = tmp1(i)
        NEXT i

        FOR i = 0 TO n - 1                    'stabilize with trapezoidal rule
          fhilf(row * n + i) = .5# * (fhilf(9 * n + i) + fhilf(8 * n + i) + h * fhilf(11 * n + i))
        NEXT i

        aufrufe = aufrufe + bufol(row) + 2

      IF row = 0 THEN GOTO 20

      'Extrapolation
      extrapol row, fhilf(), y(), n, epsabs, epsrel, x, h, h0, ahead, index, extmax, bufol()

      IF aufrufe >= fmax THEN                        ' too many calls ?
                                                     ' => stop
        fhilf(9 * n) = x
        'copyvector(fhilf(8), y, n)
        FOR i = 0 TO n - 1
          fhilf(8 * n + i) = y(i)
        NEXT i
        x = x0
        'copyvector(y, fhilf(index - 2), n)
        FOR i = 0 TO n - 1
          y(i) = fhilf((index - 2) * n + i)
        NEXT i
        ibulstoe = 1
        EXIT FUNCTION

      END IF
    END IF 'neu <> 0

    IF ahead = 0 OR neu = 0 THEN           'do we need to repeat step ?
    
      neu = 1                                  'set flag for a new call
      IF row >= Min(7, extmax - 1) THEN

      ' store differently since the extrapolation scheme has at most}
      ' 11 rows, but we allow only 8 extrapolations.

        FOR j = 0 TO 6
          'copyvector(fhilf(j), fhilf(j+1), n)
          FOR i = 0 TO n - 1
            fhilf(j * n + i) = fhilf((j + 1) * n + i)
          NEXT i
        NEXT j
      END IF
      ' accuracy could not be reached from the whole extrapolation
      ' table; repeat step with smaller step size.

      IF row >= extmax + 2 THEN
        hilf = 1# * (row - extmax + 1)
        h0 = .9# * h * (.6 ^ hilf)
        IF ABS(h0) <= 4# * MACHEPS * ABS(x0) THEN   'step too small
          ibulstoe = 2
          EXIT FUNCTION
        END IF
        row = -1
      END IF
    END IF
  GOTO 15 'end of main loop

END FUNCTION 'IBULSTOE

'Utility functions
FUNCTION Min (ia, ib)
  IF ia < ib THEN
    Min = ia
  ELSE
    Min = ib
  END IF
END FUNCTION

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