'-----------------------------------------------------------------------
' The following variables are shared among different procedures
      Public alpha, beta, gamma, delta


'-----------------------------------------------------------------------
      Sub main()
'-----------------------------------------------------------------------
' Solve the dynamics of predator-prey described by:
'   dxdt(1) =  alpha*x(1) - beta *x(1)*x(2)   ... prey
'   dxdt(2) = -gamma*x(2) + delta*x(1)*x(2)   ... predator
' Model constants: alpha, beta, gamma, delta
' Variables:
'   t = time
'   n = number of equations = 2
'-----------------------------------------------------------------------
' Program Note:
'   Call RK4 of "Numerical Recipes"
' Instructor: Nam Sun Wang
'-----------------------------------------------------------------------
      Const n = 2
      Dim x(n), dxdt(n)

' Model parameters -----------------------------------------------------
      alpha = Cells(7, 2)
      beta = Cells(8, 2)
      gamma = Cells(9, 2)
      delta = Cells(10, 2)

' Assign initial values for integration --------------------------------
'      x(1) = r0
'      x(2) = f0
      x(1) = Cells(11, 2)
      x(2) = Cells(12, 2)

' Integration limits ---------------------------------------------------
      steph = 0.1
      t = 0
      tmax = 10
      istart = 1
      iend = (tmax - t) / steph

' Save the initial conditions in the main worksheet --------------------
      irow = 15
      Cells(irow, 1) = t
      Cells(irow, 2) = x(1)
      Cells(irow, 3) = x(2)

' Calculate the dynamic values by integration --------------------------
      For i = istart To iend
'       RK4 requires the evaluation of the initial dxdt ! -- Do not forget this step
          Call DERIVS(t, x(), dxdt())
'       Actual integration with routine from "Numerical Recipes"
          Call RK4(x(), dxdt(), n, t, steph, x(), DUM)
'       End of integration step
          t = t + steph
'       Save results in the main worksheet
        Cells(irow + i, 1) = t
        Cells(irow + i, 2) = x(1)
        Cells(irow + i, 3) = x(2)
      Next i

' Plot ------------------------------------------------------------------
      ActiveSheet.ChartObjects.Add(145, 75, 300, 200).Select
      ActiveChart.ChartWizard _
        Source:=Range(Cells(irow - 1, 1), Cells(irow + 100, 3)), _
        Gallery:=xlXYScatter, Format:=6, _
        PlotBy:=xlColumns, _
        CategoryLabels:=True, SeriesLabels:=True, HasLegend:=True, _
        Title:="Predator-Prey Dynamics", _
        CategoryTitle:="Time", _
        ValueTitle:="Predator, Prey"
      ActiveChart.Axes(xlValue).TickLabels.NumberFormat = "0"
      ActiveChart.Axes(xlCategory).TickLabels.NumberFormat = "0"

      End Sub


'-----------------------------------------------------------------------
      Sub DERIVS(t, x(), dxdt())
'-----------------------------------------------------------------------
' Predator-prey dynamics
' Global variables ...
'   x(1) ... prey
'   x(2) ... predator
'   dxdt(n) ... time derivative of x variables
'   t ... time
'-----------------------------------------------------------------------

      dxdt(1) = alpha * x(1) - beta * x(1) * x(2)
      dxdt(2) = -gamma * x(2) + delta * x(1) * x(2)

      End Sub


'-----------------------------------------------------------------------
' Stuff directly from Numerical Recipes
'-----------------------------------------------------------------------
Sub RK4(Y(), DYDX(), n, x, H, YOUT(), DUM)
ReDim YT(n), DYT(n), DYM(n)
HH = H * 0.5
H6 = H / 6!
XH = x + HH
For i = 1 To n
  YT(i) = Y(i) + HH * DYDX(i)
Next i
Call DERIVS(XH, YT(), DYT())
For i = 1 To n
  YT(i) = Y(i) + HH * DYT(i)
Next i
Call DERIVS(XH, YT(), DYM())
For i = 1 To n
  YT(i) = Y(i) + H * DYM(i)
  DYM(i) = DYT(i) + DYM(i)
Next i
Call DERIVS(x + H, YT(), DYT())
For i = 1 To n
  YOUT(i) = Y(i) + H6 * (DYDX(i) + DYT(i) + 2! * DYM(i))
Next i
Erase DYM, DYT, YT
End Sub