' these are the initial values used in the primary
' example of section 8.6 of Boyce & DiPrima'a
' "Elementary Differential Equations and Boundary Value Problems"

' we have an ODE of the form
'     y' = 1 - x + 4 * y
' we want to solve for values of y over the interval from
' x = 0 to x = 1.  When x = 0, y = 1.
' We are going to break this up into 100 subintervals.

' first, we get a RungeKutta4 instance:

set rk4 = New RungeKutta4

' it doesn't matter in what order we set the properties;
' we just need to get them all set before we solve.

' first we'll tell it what equation we are solving:
rk4.ODE = "1 - x + 4 * y"

' next, we enter the initial conditions (x=0, y=1).  Note that
' our "x" is "t" in the RungeKutta class:
rk4.t0 = 0: rk4.y0 = 1

' ..and we will be solving for y as x goes to 1:

rk4.tmax = 1

' Finally we want to test over 100 intervals, so:
rk4.intervals = 100

' solving the equation returns an array of solution points in s
s = rk4.solve

' we can join each one for easier display
for i = 0 to UBound(s)
s(i) = Join(s(i), vbTab)
next

' then do another join to allow single-call output of the table
s = Join(s, vbLf)
WScript.Echo s

Class RungeKutta4

' a VBScript Implementation of a 4th-order Runge-Kutta solution to an
' ordinary differential equation
' based on reference from Section 8.6 of Boyce & DiPrima's 4th Edition
' "Elementary Differential Equations and Boundary Value Problems"
'
' Basic Runge-Kutta implementation
' provide the following values for the calculation:
' + initial values (t0, y0)
' + final t value
' + number of steps between t0 and t
' + ODE - the first order differential to be solved.
'
' the calculate method will return an array with an element
' for each step from t0 to t; each element will consist of
' the calculated t and y for the step.
'
' CAVEATS ON USE:
' All numerical methods have limitations.  It is best to ensure continuity
' of the range over the solution domain when using them; also, since
' numerical methods will always exhibit rounding errors, at a certain
' point using smaller step sizes to decrease discretization errors will
' cause even greater roundoff discussion of the issues involved with this.

dim m_t0, m_y0, m_h, m_tmax, m_intervals, m_ymax

Property Let t0(val)
m_t0 = val
End Property

Property Let y0(val)
m_y0 = val
End Property

Property get ymax
' found in calculate
' this is the value of y at the end of the interval of interest;
' useful if you don't need intermediary results on the interval.
ymax = m_ymax
end property

Property Let ODE(equation)
dim fn
fn = Array("function RHS(t,y)", "RHS = CDbl( " & equation & " )",  _
"End function")
ExecuteGlobal Join(fn, vbLf)
end property

Property Let tmax(val)
m_tmax = val
End Property

property let intervals(num)
m_intervals = num
m_h = (m_tmax - m_t0)/num
end property

function solve
dim rk(), i, kk(4)
redim rk(m_intervals)
t = m_t0
y = m_y0
rk(0) = Array(m_t0, m_y0)
for i = 1 to m_intervals
kk(1) = RHS(t, y)
kk(2) = RHS(t + .5*m_h, y + .5*m_h*kk(1))
kk(3) = RHS(t + .5*m_h, y + .5*m_h*kk(2))
kk(4) = RHS(t + m_h, y + m_h*kk(3))
y = y + m_h/6 * ( kk(1) + 2*kk(2) + 2*kk(3) + kk(4) )
t = m_t0 + (m_tmax - m_t0)/m_intervals*i
rk(i) = Array(t,y)
next
m_ymax = y
solve = rk
end function

End Class