Attribute VB_Name = "RKF45"

' Copyright (c) 1997 C. David Eagle Jr.

DefInt I-N
DefDbl A-H, O-Z

Static Sub RKF45(NEQ, TI, TF, H, TETOL, Y())

    ' Runge-Kutta-Fehlberg 4(5) subroutine

    ' Fourth-order solution with fifth-order error control
 
    ' Input

    '  NEQ   = number of equations in system
    '  TI    = initial simulation time
    '  TF    = final simulation time
    '  H     = guess for step size
    '  TETOL = truncation error tolerance
    '  Y()   = initial integration vector ( NEQ rows )

    ' Output

    '  Y()  = final integration vector ( NEQ rows )

    ' NOTE: requires SUB DERIVATIVE

    If (ICOEFF = 0) Then

       ' dimension and define integration coefficients
      
       ReDim ALPH(6), CH(6), CE(6), BETA(6, 5)

       ALPH(2) = 1# / 4#
       ALPH(3) = 3# / 8#
       ALPH(4) = 12# / 13#
       ALPH(5) = 1#
       ALPH(6) = 0.5
      
       BETA(2, 1) = 1# / 4#
       BETA(3, 1) = 3# / 32#
       BETA(3, 2) = 9# / 32#
       BETA(4, 1) = 1932# / 2197#
       BETA(4, 2) = -7200# / 2197#
       BETA(4, 3) = 7296# / 2197#
       BETA(5, 1) = 439# / 216#
       BETA(5, 2) = -8#
       BETA(5, 3) = 3680# / 513#
       BETA(5, 4) = -845# / 4104#
       BETA(6, 1) = -8# / 27#
       BETA(6, 2) = 2#
       BETA(6, 3) = -3544# / 2565#
       BETA(6, 4) = 1859# / 4104#
       BETA(6, 5) = -11# / 40#

       CH(1) = 16# / 135#
       CH(2) = 0#
       CH(3) = 6656# / 12825#
       CH(4) = 28561# / 56430#
       CH(5) = -9# / 50#
       CH(6) = 2# / 55#

       CE(1) = -1# / 360#
       CE(2) = 0#
       CE(3) = 128# / 4275#
       CE(4) = 2197# / 75240#
       CE(5) = -1# / 50#
       CE(6) = -2# / 55#

       ICOEF = 1
    End If

    ReDim YDOT(NEQ), YWRK(NEQ), F(NEQ, 6)

    ' compute integration "direction"

    SDT = Sgn(TF - TI)
    DT = Abs(H) * SDT
    DT1 = DT
    EDS = TETOL / 20#

    Do
       ' load "working" time and integration vector

       TWRK = TI
      
       For I = 1 To NEQ
           YWRK(I) = Y(I)
       Next I

       ' evaluate system of differential equations

       Call DERIVATIVE(TI, Y(), YDOT())

       For I = 1 To NEQ
           F(I, 1) = YDOT(I)
       Next I

       If (Abs(DT) < 0.0000000001) Then Exit Do

       ' compute solution

       For K = 2 To 6
           KK = K - 1
           For I = 1 To NEQ
               TEMP = 0#
               For J = 1 To KK
                   TEMP = TEMP + BETA(K, J) * F(I, J)
               Next J
               Y(I) = YWRK(I) + DT * TEMP
           Next I
           TI = TWRK + ALPH(K) * DT
           Call DERIVATIVE(TI, Y(), YDOT())
           For J = 1 To NEQ
               F(J, K) = YDOT(J)
           Next J
       Next K
      
       For I = 1 To NEQ
           TEMP = 0#
           For L = 1 To 6
               TEMP = TEMP + CH(L) * F(I, L)
           Next L
           Y(I) = YWRK(I) + DT * TEMP
       Next I

       ' truncation error calculations

       TER = 0#

       For I = 1 To NEQ
           A = Abs(Y(I))
          
           If (A < TETOL) Then A = TETOL

           TEI = Abs(CE(1) * F(I, 1) + CE(3) * F(I, 3) + CE(4) * F(I, 4) + CE(5) * F(I, 5) + CE(6) * F(I, 6)) / A
          
           If (TEI > TER) Then TER = TEI
       Next I
       
       TER = TER * Abs(DT) + 1E-20

       DT1 = DT

       Q = EDS / TER

       ' compute new step size

       DT = DT1 * Q ^ (1# / 5#)

       If (TER > TETOL) Then
          ' reject current step
          TI = TWRK
          For I = 1 To NEQ
              Y(I) = YWRK(I)
          Next I
       Else
          ' accept current step
          TI = TWRK + DT1
          If (Abs(DT) > Abs(TF - TI)) Then DT = TF - TI
       End If
      
    Loop

    ' erase working arrays

    Erase F, YDOT, YWRK

End Sub