'***********************************************************************
'* Integrate a System of Ordinary Differential Equations By the *
'* Runge-Kutta-Fehlberg method (double precision) *
'* ------------------------------------------------------------------- *
'* REFERENCE: H A Watts and L F Shampine, *
'* Sandia Laboratories, *
'* Albuquerque, New Mexico. *
'* ------------------------------------------------------------------- *
'* SAMPLE RUN: *
'* *
'* PROGRAM TRKF45 *
'* Demonstrate the RKF45 ODE integrator. *
'* *
'* TEST01 *
'* Solve a scalar equation: *
'* *
'* Y' = 0.25 * Y * ( 1 - Y / 20 ) *
'* *
'* T Y Y_Exact Error *
'* *
'* 0.00000 1.00000 1.00000 0.0000000 *
'* 4.00000 2.50321 2.50322 -0.0000087 *
'* 8.00000 5.60007 5.60009 -0.0000193 *
'* 12.00000 10.27774 10.27773 -0.0000069 *
'* 16.00000 14.83682 14.83683 -0.0000038 *
'* 20.00000 17.73017 17.73017 -0.0000084 *
'* *
'* TEST02 *
'* Solve a vector equation: *
'* *
'* Y'(1) = Y(2) *
'* Y'(2) = - Y(1) *
'* *
'* T Y1 Y2 *
'* *
'* 0.00000 1.00000 0.00000 *
'* 0.52360 0.86603 -0.50000 *
'* 1.04720 0.50000 -0.86603 *
'* 1.57080 0.00000 -1.00000 *
'* 2.09440 -0.50000 -0.86603 *
'* 2.61799 -0.86603 -0.50000 *
'* 3.14159 -1.00000 -0.00000 *
'* 3.66519 -0.86603 0.50000 *
'* 4.18879 -0.50000 0.86603 *
'* 4.71239 -0.00000 1.00001 *
'* 5.23599 0.50000 0.86604 *
'* 5.75959 0.86604 0.50001 *
'* 6.28319 1.00002 0.00000 *
'* *
'* TEST03 *
'* Solve a vector equation: *
'* *
'* Y'(1) = Y(2) *
'* Y'(2) = Y(3) *
'* Y'(3) = Y(4) *
'* Y'(4) = Y(5) *
'* Y'(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4)^3) / (9 * Y(3)^2) *
'* *
'* T Y1 Y2 Y3 Y4 Y5 *
'* *
'* 0.00000 1.00000 1.00000 1.00000 1.00000 1.00000 *
'* 0.13636 1.14610 1.14609 1.14587 1.14068 1.05604 *
'* 0.27273 1.31354 1.31340 1.31128 1.28532 1.05248 *
'* 0.40909 1.50538 1.50460 1.49612 1.42333 0.95209 *
'* 0.54545 1.72508 1.72223 1.69844 1.53859 0.71111 *
'* 0.68182 1.97638 1.96840 1.91370 1.60897 0.28809 *
'* 0.81818 2.26328 2.24438 2.13400 1.60781 -0.33918 *
'* 0.95455 2.58984 2.55011 2.34770 1.50801 -1.15027 *
'* 1.09091 2.96003 2.88369 2.53985 1.28946 -2.06094 *
'* 1.22727 3.37739 3.24105 2.69376 0.94820 -2.92046 *
'* 1.36364 3.84475 3.61589 2.79372 0.50379 -3.54302 *
'* 1.50000 4.36396 4.00000 2.82843 -0.00000 -3.77124 *
'* ------------------------------------------------------------------- *
'* *
'* Basic Release 1.1 By J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'* *
'* Release 1.1: added test #3. *
'***********************************************************************
' LIST OF SUBROUTINES:
' 400: User defined calculation of derivatives (examples #1 and #2)
' 500: Calculate exact solution (example #1)
' 600: Define example #1
' 700: Define example #2
' 800: Define example #3
' 1000: Fehlberg subroutine (one step)
' 1200: Emulation of Fortran routine Sign (with integers)
' 1210: Emulation of Fortran routine Sign (with real numbers)
' 1300: Emulation of Fortran routine Max (with real numbers)
' 1310: Emulation of Fortran routine Min (with real numbers)
' 2000: Runge-Kutta-Fehlberg method (double precision)
'------------------------------------------------------------
DefDbl A-H, O-Z
DefInt I-N
NEQ = 5 'Maximum number of equations
PI = 4# * Atn(1#) 'Constant PI
EPSILON = 2.22E-16 'Small number
Dim y(NEQ), yp(NEQ)
'auxiliary variables used by 400
Dim yy(NEQ), yyp(NEQ)
'work space used by successive calls to 1000 and 2000
Dim f1(NEQ), f2(NEQ), f3(NEQ), f4(NEQ), f5(NEQ)
Cls
Print
Print " PROGRAM TRKF45"
Print " Demonstrate the RKF45 ODE integrator."
GoSub 600 'call test01
Print
INPUT "Any key to continue... ", R$
GoSub 700 'call test02
Print
INPUT "Any key to continue... ", R$
GoSub 800 'call test03
Print
INPUT "Any key to continue... ", R$
End 'of main program
'User defined system of diff. equations
400 'Subroutine f(tt, yy, yyp)
'--------------------------------------------------------------------
'
' F evaluates the derivative for the ODE (TESTS #1 and #2).
'
'-------------------------------------------------------------------}
If num = 1 Then
yyp(1) = 0.25 * yy(1) * (1# - yy(1) / 20#)
ElseIf num = 2 Then
yyp(1) = yy(2): yyp(2) = -yy(1)
Else
yyp(1) = yy(2)
yyp(2) = yy(3)
yyp(3) = yy(4)
yyp(4) = yy(5)
yyp(5) = (45# * yy(3) * yy(4) * yy(5) - 40# * yy(4) ^ 3) / (9# * yy(3) ^ 2)
End If
Return
500 'function yexact(tt)
'--------------------------------------------------------------------
'
' YEXACT evaluates the exact solution of the ODE (For TEST #1).
'
'--------------------------------------------------------------------
yexact = 20# / (1# + 19# * Exp(-0.25 * tt))
Return
600 'Subroutine test01
'--------------------------------------------------------------------
'
' TEST01 solves a scalar ODE in double precision.
'
'--------------------------------------------------------------------
num = 1 'example #1
neqn = 1 'one equation
F$ = "####.#####"
G$ = "####.#######"
Print
Print " TEST01"
Print " Solve a scalar equation:"
Print
Print " Y' = 0.25 * Y * ( 1 - Y / 20 )"
Print
abserr = 0.000001
relerr = 0.000001
iflag = 1
tstart = 0#
tstop = 20#
nstep = 5
tout = 0#
y(1) = 1#
Print " T Y Y_Exact Error"
Print
tt = tout: GoSub 500 'calculate yexact(tout)
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; yexact;
Print " ";: Print USING; G$; y(1) - yexact
For istep = 1 To nstep
t = ((nstep - istep + 1) * tstart + (istep - 1) * tstop) / nstep
tout = ((nstep - istep) * tstart + (istep) * tstop) / nstep
GoSub 2000 'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,
'nfe,kop,init,jflag,kflag)
tt = tout: GoSub 500 'calculate yexact(tout)
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; yexact;
Print " ";: Print USING; G$; y(1) - yexact
Next istep
Return
700 'Subroutine test02
'--------------------------------------------------------------------
'
' TEST02 solves a vector ODE.
'
'--------------------------------------------------------------------
num = 2 'Example #2
neqn = 2 '2 equations
Print
Print " TEST02"
Print " Solve a vector equation:"
Print
Print " Y'(1) = Y(2)"
Print " Y'(2) = - Y(1)"
abserr = 0.000001
relerr = 0.000001
iflag = 1
tstart = 0#
tstop = 2# * PI
nstep = 12
tout = 0#
y(1) = 1#
y(2) = 0#
Print
Print " T Y1 Y2"
Print
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; y(2)
For istep = 1 To nstep
t = ((nstep - istep + 1) * tstart + (istep - 1) * tstop) / nstep
tout = ((nstep - istep) * tstart + (istep * tstop) / nstep)
GoSub 2000 'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,
'nfe,kop,init,jflag,kflag)
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; y(2)
Next istep
Return
800 'Subroutine test03
'--------------------------------------------------------------------
'
' TEST02 solves a vector ODE.
'
'--------------------------------------------------------------------
num = 3 'Example #3
neqn = 5 '5 equations
Print
Print " TEST03"
Print " Solve a vector equation:"
Print
Print " Y'(1) = Y(2)"
Print " Y'(2) = Y(3)"
Print " Y'(3) = Y(4)"
Print " Y'(4) = Y(5)"
Print " Y'(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4)^3) / (9 * Y(3)^2)"
abserr = 0.000001
relerr = 0.000001
iflag = 1
tstart = 0#
tstop = 1.5
nstep = 11
tout = 0#
y(1) = 1#
y(2) = 1#
y(3) = 1#
y(4) = 1#
y(5) = 1#
Print
Print " T Y1 Y2 Y3 Y4 Y5"
Print
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; y(2);
Print " ";: Print USING; F$; y(3);
Print " ";: Print USING; F$; y(4);
Print " ";: Print USING; F$; y(5)
For istep = 1 To nstep
t = ((nstep - istep + 1) * tstart + (istep - 1) * tstop) / nstep
tout = ((nstep - istep) * tstart + (istep * tstop) / nstep)
GoSub 2000 'call rkfs (neqn,y,t,tout,relerr,abserr,iflag,yp,h,f1..f5,savr,save,
'nfe,kop,init,jflag,kflag)
Print " ";: Print USING; F$; tout;
Print " ";: Print USING; F$; y(1);
Print " ";: Print USING; F$; y(2);
Print " ";: Print USING; F$; y(3);
Print " ";: Print USING; F$; y(4);
Print " ";: Print USING; F$; y(5)
Next istep
Return
'************************************************************************
'* Collection of Basic subroutines to Integrate a System of Ordinary *
'* Differential Equations By the Runge-Kutta-Fehlberg method (in double *
'* precision). *
'************************************************************************
1000 'Subroutine fehl(neqn, y, t, h, yp, f1, f2, f3, f4, f5, s)
'************************************************************************
'
' FEHL takes one Fehlberg fourth-fifth order step (double precision).
'
' Discussion:
'
' FEHL integrates a system of NEQN first order ordinary differential
' equations of the form
' dY(i)/dT = F(T,Y(1),---,Y(NEQN))
' where the initial values Y and the initial derivatives
' YP are specified at the starting point T.
'
' FEHL advances the solution over the fixed step H and returns
' the fifth order (sixth order accurate locally) solution
' approximation at T+H in array S.
'
' The formulas have been grouped to control loss of significance.
' FEHL should be called with an H not smaller than 13 units of
' roundoff in T so that the various independent arguments can be
' distinguished.
'
' Author:
'
' H A Watts and L F Shampine,
' Sandia Laboratories,
' Albuquerque, New Mexico.
'
' RKF45 is primarily designed to solve non-stiff and mildly stiff
' differential equations when derivative evaluations are inexpensive.
'
' Parameters:
'
' Input, external F, a subroutine of the form
' Procedure f(t:double; y:VEC; var yp:VEC);
' to evaluate the derivatives.
' YP(I) = dY(I) / dT;
'
' Input, integer NEQN, the number of equations to be integrated.
'
' Input, double Y(NEQN), the current value of the dependent variable.
'
' Input, double T, the current value of the independent variable.
'
' Input, double H, the step size to take.
'
' Input, double YP(NEQN), the current value of the derivative of the
' dependent variable.
'
' Output, double F1(NEQN), double F2(NEQN), double F3(NEQN), double F4(NEQN),
' double F5(NEQN) are arrays of dimension NEQN which are needed for
' internal storage.
'
' Output, double S(NEQN), the computed estimate of the solution at T+H.
'******************************************************************************}
ch = h / 4#
For i = 1 To neqn
f5(i) = y(i) + ch * yp(i)
Next i
tt = t + ch:
For i = 1 To neqn
yy(i) = f5(i)
Next i
GoSub 400 'call f(t+ch, f5, f1)
For i = 1 To neqn
f1(i) = yyp(i)
Next i
ch = 3# * h / 32#
For i = 1 To neqn
f5(i) = y(i) + ch * (yp(i) + 3# * f1(i))
Next i
tt = t + 3# * h / 8#
For i = 1 To neqn
yy(i) = f5(i)
Next i
GoSub 400 'call f(t+3#*h/8#, f5, f2)
For i = 1 To neqn
f2(i) = yyp(i)
Next i
ch = h / 2197#
For i = 1 To neqn
f5(i) = y(i) + ch * (1932# * yp(i) + (7296# * f2(i) - 7200# * f1(i)))
Next i
tt = t + 12# * h / 13#
For i = 1 To neqn
yy(i) = f5(i)
Next i
GoSub 400 'call f(t+12#*h/13#, f5, f3)
For i = 1 To neqn
f3(i) = yyp(i)
Next i
ch = h / 4104#
For i = 1 To neqn
f5(i) = y(i) + ch * ((8341# * yp(i) - 845# * f3(i)) + (29440# * f2(i) - 32832# * f1(i)))
Next i
tt = t + h
For i = 1 To neqn
yy(i) = f5(i)
Next i
GoSub 400: 'call f(t+h, f5, f4)
For i = 1 To neqn
f4(i) = yyp(i)
Next i
ch = h / 20520#
For i = 1 To neqn
tmp = (41040# * f1(i) - 28352# * f2(i))
f1(i) = y(i) + ch * ((-6080# * yp(i) + (9295# * f3(i) - 5643# * f4(i))) + tmp)
Next i
tt = t + h / 2#
For i = 1 To neqn
yy(i) = f1(i)
Next i
GoSub 400 'call f(t+h/2#, f1, f5)
For i = 1 To neqn
f5(i) = yyp(i)
Next i
' Ready to compute the approximate solution at T+H.
ch = h / 7618050#
For i = 1 To neqn
tmp = (3953664# * f2(i) + 277020# * f5(i))
f1(i) = y(i) + ch * ((902880# * yp(i) + (3855735# * f3(i) - 1371249# * f4(i))) + tmp)
Next i
Return
'Emulation of Fortran Intrisic Functions
1200 'Function ISign(ia,ib)
If ib < 0 Then
ISign = -Abs(ia)
Else
ISign = Abs(ia)
End If
Return
1210 'Function Sign(a,b)
If B < 0 Then
Sign = -Abs(a)
Else
Sign = Abs(a)
End If
Return
1300 'Function XMax(a,b)
If a >= B Then
XMax = a
Else
XMax = B
End If
Return
1310 'Function XMin(a,b)
If a <= B Then
XMIN = a
Else
XMIN = B
End If
Return
2000 'Subroutine rkfs(neqn, y, t, tout, relerr, abserr, iflag, yp, h, f1, f2, f3, f4,f5,
' savr, save, nfe, kop, init, jflag, kflag)
'***************************************************************************************
'
' RKFS implements the Runge-Kutta-Fehlberg method (double precision).
'
' Discussion:
'
' RKFS integrates a system of first order ordinary differential
' equations as described in the comments for RKF45.
'
' The arrays yp, f1, f2, f3, f4, and f5 (of dimension neqn) and
' the variables h, savre, savae, nfe, kop, init, jflag and kflag are used
' internally by the code and appear in the call list to eliminate
' local retention of variables between calls. Accordingly, they
' should not be altered. Items of possible interest are
'
' YP - the derivative of the solution vector at T;
' H - an appropriate stepsize to be used for the next step;
' NFE - the number of derivative function evaluations.
'
' The expense is controlled by restricting the number
' of function evaluations to be approximately MAXNFE.
' As set, this corresponds to about 500 steps.
'
' REMIN is the minimum acceptable value of RELERR. Attempts
' to obtain higher accuracy with this subroutine are usually
' very expensive and often unsuccessful.
'***************************************************************************************
'Labels: 25,40,45,50,60,65,80,100,200,260
remin = 0.000000000001
maxnfe = 3000
' ihfaild, ioutput: boolean (0 or 1)
' Check the input parameters
eps = EPSILON
If neqn < 1 Then
iflag = 8
Return
End If
If relerr < 0# Then
iflag = 8
Return
End If
If abserr < 0# Then
iflag = 8
Return
End If
mflag = Abs(iflag)
If Abs(iflag) < 1 Or Abs(iflag) > 8 Then
iflag = 8
Return
End If
' Is this the first call?
If mflag = 1 Then GoTo 50
' Check continuation possibilities
If t = tout And kflag <> 3 Then
iflag = 8
Return
End If
If mflag <> 2 Then GoTo 25
' iflag = +2 or -2
If kflag = 3 Then GoTo 45
If init = 0 Then GoTo 45
If kflag = 4 Then GoTo 40
If kflag = 5 And abserr = 0# Then End
If kflag = 6 And relerr <= savr And abserr <= save Then End
GoTo 50
' iflag = 3,4,5,6,7 or 8
25 If iflag = 3 Then GoTo 45
If iflag = 4 Then GoTo 40
If iflag = 5 And abserr > 0# Then GoTo 45
' Integration cannot be continued since user did not respond to
' the instructions pertaining to iflag=5,6,7 or 8.
End
' Reset function evaluation counter
40 nfe = 0
If mflag = 2 Then GoTo 50
' Reset flag value from previous call
45 iflag = jflag
If kflag = 3 Then mflag = Abs(iflag)
' Save input iflag and set continuation flag for subsequent input checking
50 jflag = iflag
kflag = 0
' Save relerr and abserr for checking input on subsequent calls
savr = relerr
save = abserr
' Restrict relative error tolerance to be at least as large as
' 2*eps+remin to avoid limiting precision difficulties arising
' from impossible accuracy requests.
rer = 2# * EPSILON + remin
' The relative error tolerance is too small
If relerr < rer Then
relerr = rer
iflag = 3
kflag = 3
Return
End If
dt = tout - t
If mflag = 1 Then GoTo 60
If init = 0 Then GoTo 65
GoTo 80
' Initialization:
' set initialization completion indicator, init
' set indicator for too many output points, kop
' evaluate initial derivatives
' set counter for function evaluations, nfe
' evaluate initial derivatives
' set counter for function evaluations, nfe
' estimate starting stepsize.
60 init = 0
kop = 0
a = t
For i = 1 To neqn
yy(i) = y(i)
Next i
GoSub 400 'call f(a, y, yp)
For i = 1 To neqn
yp(i) = yyp(i)
Next i
nfe = 1
If t = tout Then
iflag = 2
Return
End If
65 init = 1
h = Abs(dt)
toln = 0#
For k = 1 To neqn
tol = relerr * Abs(y(k)) + abserr
If tol > 0# Then
toln = tol
ypk = Abs(yp(k))
If ypk * h ^ 5 > tol Then
h = (tol / ypk) ^ 0.2
End If
End If
Next k
If toln <= 0# Then h = 0#
'h = XMax(h, 26# * eps * XMax(abs(t), abs(dt)))
If Abs(t) > Abs(dt) Then
tmp = Abs(t)
Else
tmp = Abs(dt)
End If
If h < 26# * eps * tmp Then
h = 26# * eps * tmp
End If
ia = 2: ib = iflag: GoSub 1200
jflag = ISign
' Set stepsize for integration in the direction from T to TOUT
80 a = h: B = dt: GoSub 1210
h = Sign
' Test to see if RKF45 is being severely impacted by too many output points.
If Abs(h) >= 2# * Abs(dt) Then kop = kop + 1
' Unnecessary frequency of output
If kop = 100 Then
kop = 0
iflag = 7
Return
End If
' If too close to output point, extrapolate and return
If Abs(dt) <= 26# * eps * Abs(t) Then
For i = 1 To neqn
y(i) = y(i) + dt * yp(i)
Next i
a = tout
For i = 1 To neqn
yy(i) = y(i)
Next i
GoSub 400 'call f(a, y, yp)
For i = 1 To neqn
yp(i) = yyp(i)
Next i
nfe = nfe + 1
t = tout
iflag = 2
Return
End If
' Initialize output point indicator
ioutput = 0
' To avoid premature underflow in the error tolerance function,
' scale the error tolerances.
scale = 2# / relerr
ae = scale * abserr
' Step by step integration
100: ihfaild = 0
' Set smallest allowable stepsize
hmin = 26# * eps * Abs(t)
' Adjust stepsize if necessary to hit the output point.
' Look ahead two steps to avoid drastic changes in the stepsize and
' thus lessen the impact of output points on the code.
dt = tout - t
If Abs(dt) >= 2# * Abs(h) Then GoTo 200
' The next successful step will complete the integration to the output point.
If Abs(dt) <= Abs(h) Then
ioutput = 1
h = dt
GoTo 200
End If
h = 0.5 * dt 'reduce step
' Core integrator for taking a single step
'
' The tolerances have been scaled to avoid premature underflow in
' computing the error tolerance function ET.
' To avoid problems with zero crossings, relative error is measured
' using the average of the magnitudes of the solution at the
' beginning and end of a step.
' The error estimate formula has been grouped to control loss of
' significance.
'
' To distinguish the various arguments, H is not permitted
' to become smaller than 26 units of roundoff in T.
' Practical limits on the change in the stepsize are enforced to
' smooth the stepsize selection process and to avoid excessive
' chattering on problems having discontinuities.
' To prevent unnecessary failures, the code uses 9/10 the stepsize
' it estimates will succeed.
'
' After a step failure, the stepsize is not allowed to increase for
' the next attempted step. This makes the code more efficient on
' problems having discontinuities and more effective in general
' since local extrapolation is being used and extra caution seems
' warranted.
'
' Test number of derivative function evaluations.
' If okay, try to advance the integration from T to T+H.
' Too many function calls.
200 If nfe > maxnfe Then
iflag = 4
kflag = 4
Return
End If
' Advance an approximate solution over one step of length H.
GoSub 1000 'call fehl(neqn, y, t, h, yp, f1, f2, f3, f4, f5, f1)
nfe = nfe + 5
' Compute and test allowable tolerances versus local error estimates
' and remove scaling of tolerances. Note that relative error is
' measured with respect to the average of the magnitudes of the
' solution at the beginning and end of the step.
eeoet = 0#
For k = 1 To neqn
et = Abs(y(k)) + Abs(f1(k)) + ae
If et <= 0# Then
iflag = 5
Return
End If
tmp = (22528# * f2(k) - 27360# * f5(k))
ee = Abs((-2090# * yp(k) + (21970# * f3(k) - 15048# * f4(k))) + tmp)
'eeoet = XMax(eeoet, ee/et)
a = eecet: B = ee / et: GoSub 1300: eecet = XMax
Next k
esttol = ABS(h) * eeoet * scale / 752400#
If esttol <= 1# Then GoTo 260
' Unsuccessful step. Reduce the stepsize, try again.
' The decrease is limited to a factor of 1/10.
ihfaild = 1
ioutput = 0
If esttol < 59049# Then
s = 0.9 / (esttol ^ 0.2)
Else
s = 0.1
End If
h = s * h
If Abs(h) < hmin Then
iflag = 6
kflag = 6
Return
Else
GoTo 200
End If
' Successful step. Store solution at T+H and evaluate derivatives there.
260 t = t + h
For i = 1 To neqn
y(i) = f1(i)
Next i
a = t
For i = 1 To neqn
yy(i) = y(i)
Next i
GoSub 400 'call f(a, y, yp)
For i = 1 To neqn
yp(i) = yyp(i)
Next i
nfe = nfe + 1
' Choose next stepsize. The increase is limited to a factor of 5.
' If step failure has just occurred, next stepsize is not allowed to increase.
If esttol > 0.0001889568 Then
s = 0.9 / esttol ^ 0.2
Else
s = 5#
End If
If ihfaild <> 0 Then
's = XMin(s, 1#)
a = s: B = 1#: GoSub 1310: s = XMIN
End If
'h = Sign(XMax(s * abs(h), hmin), h)
a = s * Abs(h): B = hmin: GoSub 1300
a = XMax: B = h: GoSub 1210: h = Sign
' End of core integrator
' Should we take another step?
If ioutput <> 0 Then
t = tout
iflag = 2
End If
If iflag > 0 Then GoTo 100
' Integration successfully completed
' ne-step mode
iflag = -2
Return
'end of file trkf45.bas
