DECLARE SUB odeint (ystart#(), nvar%, t1#, t2#, eps#, h1#, hmin#, nok%, nbad%)
DECLARE SUB Derivs (t#, XT#(), XF#())
DECLARE SUB rk4 (y#(), dydt#(), n%, t#, h#, yout#())
DECLARE SUB rkqc (y#(), dydt#(), n%, t#, htry#, eps#, yscal#(), hdid#, hnext#)
'************************************************************
'*            TEST PROGRAM OF SUBROUTINE ODEINT             *
'*       (Runge-Kutta Method with Time Step Control)        *
'* -------------------------------------------------------- *
'* SAMPLE RUN:                                              *
'* Integrate ODE First Order System:                        *
'*   Y1' = Y2                                               *
'*   Y2' = (3-C)*Y1+2*Y4-q1                                 *
'*   Y3' = Y4                                               *
'*   Y4' = -2*Y2-C*Y3                                       *
'*   Y5' = Y6                                               *
'*   Y6' = -(1+C)*Y5                                        *
'* from t1=0 to t2=10 with additional data:                 *
'*   Y start = (0,0,0,0,1,0)                                *
'*   Starting integration step = 0.01                       *
'*   Desired precision = 1E-8                               *
'*   Minimum integration step = 0.001                       *
'*   C=1, q1=2                                              *
'*                                                          *
'* At time =    10.0000000000000                            *
'* Y(1) = -0.121038E+04                                     *
'* Y(2) = -0.908048E+03                                     *
'* Y(3) =  0.116298E+04                                     *
'* Y(4) =  0.870537E+03                                     *
'* Y(5) = -0.496866E-02                                     *
'* Y(6) = -0.141420E+01                                     *
'*                                                          *
'* Final time step = 2.175370624046558E-002                 *
'*                                                          *
'*                Quick Basic Release By J-P Moreau, Paris. *
'*                            (www.jpmoreau.fr)             *
'************************************************************
'Program Tequdif
DEFDBL A-H, O-Z
DEFINT I-N

CONST MAXVAR = 10
DIM y(MAXVAR)

  nvar = 6
  FOR i = 1 TO nvar
    y(i) = 0#
  NEXT i
  y(5) = 1#
  t1 = 0#
  t2 = 15#
  h = .01#
  dt = .5#

  F$ = "###.##"
  F1$ = "#######.####"

  t = t1
  CLS
  PRINT
  PRINT "  Time      Y1          Y2          Y3          Y4          Y5         Y6"
  PRINT "------------------------------------------------------------------------------"
  PRINT USING F$; t;
  FOR i = 1 TO nvar
    PRINT USING F1$; y(i);
  NEXT i
  PRINT

  WHILE t < t2
    odeint y(), nvar, t, t + dt, .0000000001#, h, .0001#, nok, nbad
    PRINT USING F$; t + dt;
    FOR i = 1 TO nvar
      PRINT USING F1$; y(i);
    NEXT i
    PRINT
    t = t + dt
  WEND
  PRINT "------------------------------------------------------------------------------"
  PRINT
  PRINT " Final time step = "; h

END 'of main program


'first order ODE system to integrate
SUB Derivs (t, XT(), XF())
CONST C = 1#, q1 = 2#
  XF(1) = XT(2)
  XF(2) = (3# - C) * XT(1) + 2# * XT(4) - q1
  XF(3) = XT(4)
  XF(4) = -2# * XT(2) - C * XT(3)
  XF(5) = XT(6)
  XF(6) = -(1# + C) * XT(5)
END SUB

SUB odeint (ystart(), nvar, t1, t2, eps, h1, hmin, nok, nbad)
'================================================================
' From Pascal Version 1.2 - 12/31/1991 By J-P Dumont, France.
'================================================================
'  This procedure integrates the 1st order differential system:
'       dy/dt = F(y,t)
'  where y and F are vectors of size nvar, between times
'  t1 and t2.
'
'  INPUTS:
'  ystart= begin coordinates vector and speeds
'  nvar  = number of equations
'  t1    = begin integration time
'  t2    = end integration time
'          t2 may be > or < t1
'  eps   = absolute required precision for solution
'  h1    = time increment proposed at beginning (try value)
'  hmin  = minimum time increment
'
'  OUTPUTS:
'  ystart= end coordinates vector and speeds
'  nok   = number of unchanged time steps
'  nbad  = number of modified time steps
'==================================================================
' Programs using procedure RKQC must provide access to:
'
' Subroutine Derivs(t, y(), dydt())
' which returns the values of derivatives dydt at point t, knowing
' both t and the values of functions y.
'
'==================================================================

CONST maxstp = 10000, two = 2#, zero = 0#, tiny = 9.999999999999999D-21

DIM yscal(MAXVAR), y(MAXVAR), dydt(MAXVAR)

     t = t1
     IF t2 > t1 THEN
       h = ABS(h1)
     ELSE
       h = -ABS(h1)
     END IF
     nok = 0
     nbad = 0
     FOR i = 1 TO nvar
       y(i) = ystart(i)
     NEXT i
     tsav = t - dtsave * two
     FOR nstp = 1 TO maxstp
         Derivs t, y(), dydt()
         FOR i = 1 TO nvar
           yscal(i) = ABS(y(i)) + ABS(dydt(i) * h)
         NEXT i
         IF ((t + h - t2) * (t + h - t1)) > zero THEN h = t2 - t
         rkqc y(), dydt(), nvar, t, h, eps, yscal(), hdid, hnext
         IF hdid = h THEN
           nok = nok + 1
         ELSE
           nbad = nbad + 1
         END IF
         IF ((t - t2) * (t2 - t1)) >= zero THEN
            FOR i = 1 TO nvar
              ystart(i) = y(i)
            NEXT i
            GOTO 99   'it's over
         END IF
         IF ABS(hnext) < hmin THEN
           PRINT " Time step too small!"
           nok = -1   'error flag
           INPUT "", RR$
           GOTO 99
         END IF
         h = hnext
     NEXT nstp
     INPUT " Pause in Sub ODEINT - too many time steps!", RR$
99   h1 = h'return last time step
END SUB

SUB rk4 (y(), dydt(), n, t, h, yout())
'=================================================================
'  From Pascal Version 1.1 - 12/31/1991 By J-P Dumont, France.
'=================================================================
' This procedure integrates the 1st order differential system:
'      dy/dt = F(y,t)
' by a fourth order Runge-Kutta method to advance the solution
' on interval h of the independant variable t:
'-----------------------------------------------------------------
' INPUTS:
' y     = State vector at begin of integration step
' dydt  = its derivative at the same point.
' n     = number of equations of system
' t     = i.s. value at begin of time step
' h     = time step
'------------------------------------------------------------------
' OUTPUT:
' yout  = State vector at end of integration step.
'------------------------------------------------------------------
' Programs using procedure RK4 must provide access to:
'
' Subroutine Derivs(t, y(), dydt())
' which returns the values of derivatives dydt at point t, knowing
' both t and the values of functions y.
'
'------------------------------------------------------------------
DIM dym(MAXVAR), dyt(MAXVAR), yt(MAXVAR)
  hh = .5 * h
  h6 = h / 6
  th = t + hh
  FOR i = 1 TO n
    yt(i) = y(i) + hh * dydt(i)
  NEXT i
  Derivs th, yt(), dyt()
  FOR i = 1 TO n
    yt(i) = y(i) + hh * dyt(i)
  NEXT i
  Derivs th, yt(), dym()
  FOR i = 1 TO n
    yt(i) = y(i) + h * dym(i)
    dym(i) = dyt(i) + dym(i)
  NEXT i
  Derivs t + h, yt(), dyt()
  FOR i = 1 TO n
    yout(i) = y(i) + h6 * (dydt(i) + dyt(i) + 2# * dym(i))
  NEXT i
END SUB

SUB rkqc (y(), dydt(), n, t, htry, eps, yscal(), hdid, hnext)
'=====================================================================
' From Pascal Version 1.2 - 03/26/1993 By J-P Dumont, France.
'=====================================================================
' Runge-Kutta integration step with control of truncation local error
' to obtain a required precision and adjust time step consequently
'---------------------------------------------------------------------
' INPUTS:
' y      = State vector of size n
' dydt   = its derivative at begin value of independant variable, t
' n      = number of equations of system
' t      = begin value of independant variable
' htry   = time step proposed as a try
' eps    = precision requirement:
'          Max (ycalc[i] - yvrai[i])/yscal[i] < eps
' yscal  = normalization vector of solution.
'---------------------------------------------------------------------
' OUTPUTS:
'  y     = end state vector
'  t     = i.s. end value
'  hdid  = actual time step 
'  hnext = time step advised for the next integration step
'---------------------------------------------------------------------
' Programs using subroutine RKQC must provide access to:
'
' Subroutine Derivs(t, y(), dydt())
' which returns the values of derivatives dydt at point t, knowing
' both t and the values of functions y.
'
'=====================================================================

CONST pgrow = -.2, pshrnk = -.25, fcor = .06666666#'1/15
CONST un = 1#, safety = .9#, errcon = 5.999999999999999D-04, tiny = 9.999999999999999D-21

DIM dysav(MAXVAR), ysav(MAXVAR), ytemp(MAXVAR)

         tsav = t         'Save begin time
         FOR i = 1 TO n
           ysav(i) = y(i)
           dysav(i) = dydt(i)
         NEXT i
         h = htry         'define increment for a try value
1        hh = .5 * h      'take 2 half time steps
         rk4 ysav(), dysav(), n, tsav, hh, ytemp()
         t = tsav + hh
         Derivs t, ytemp(), dydt()
         rk4 ytemp(), dydt(), n, t, hh, y()
         t = tsav + h
         IF t = tsav THEN
           PRINT " Pause in Subroutine RKQC"
           PRINT " Increment too small of independant variable"
           INPUT " Press any key to continue... ", RR$
         END IF
         rk4 ysav(), dysav(), n, tsav, h, ytemp()
         errmax = 0#      'Evaluate error
         temp = 0#
         FOR i = 1 TO n
           ytemp(i) = y(i) - ytemp(i)    'ytemp = estimated error
           IF yscal(i) > tiny THEN temp = ABS(ytemp(i) / yscal(i))
           IF errmax < temp THEN errmax = temp
         NEXT i
         errmax = errmax / eps    'real error / requirement
         IF errmax > un THEN      'Error too big, reduce h
           h = safety * h * EXP(pshrnk * LOG(errmax))
           GOTO 1                 'start again
         ELSE                     'the step has been a success!
           hdid = h               'Calculate next time step
           IF errmax > errcon THEN
             hnext = safety * h * EXP(pgrow * LOG(errmax))
           ELSE
             hnext = 4# * h
           END IF
         END IF
         FOR i = 1 TO n
           y(i) = y(i) + ytemp(i) * fcor
         NEXT i
END SUB