' Dimension an array for the results
ReDim YA2(0 To M - 1, 0 To N - 1)

' Set up the "State" object to transfer data between the AlgLib code and your code to evealuate the differential equations
' See AlgLib code comments or manual for more details
Call ODESolverRKCK(YA(), N, XA2, M, Eps, Step, State)

' Loop through ODESolver and evaluation of DE until all steps have converged
Rtn = True
Do While Rtn = True
Rtn = ODESolverIteration(State)
' Either hard code evaluation of the Diff Equn(s) or call a VBA function with the name specefied in FuncName as shown below
State.DY = Application.Run(FuncName, State.X, State.Y, CoeffA)
Loop

' Extract the results array from the State function
Call ODESolverResults(State, M, XA2, YA2, Rep)
ODE = YA2

Function ODEFunc1(X As Double, Y As Variant, CoeffA As Variant) As Variant
Dim ResA(0 To 0) As Double
ResA(0) = CoeffA(1, 1) * Y(0)
ODEFunc1 = ResA
End Function

dy0/dx = f0(y0, ..., y8, x)
dy1/dx = f1(y0, ..., y8, x)
...
dy8/dx = f8(y0, ..., y8, x)