#
# Borrowed and modified from: http://www.theresearchkitchen.com/archives/679
#
runge_kutta = function(dy.dx=function(x,y){}, h=1E-7, y0=1, start=0, end=1){
    nsteps = (end-start)/h
    xs = numeric(nsteps+1)
    ys = numeric(nsteps+1)
    xs[1] = start
    ys[1] = y0
    for (i in 1:nsteps){
        x = start + (i-1)*h
        y = ys[i]
        k1 = dy.dx(x, y)
        k2 = dy.dx(x+0.5*h, y + 0.5*h*k1)
        k3 = dy.dx(x+0.5*h, y + 0.5*h*k2)
        k4 = dy.dx(x+h, y + h*k3)
        ys[i+1] = ys[i] + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
        xs[i+1] = x + h
    }
    list( xs=xs, ys=ys )
}

#
# An example of how to use use this function is:
#
# dy.dx = function(x,y) { 3*x - y + 8 }
#
# res = runge_kutta(dy.dx, start=0, end=0.5, h=0.1, y0=3)
#
# res$ys
# 
# [1] 3.000000 3.490325 3.962538 4.418363 4.859359 5.286938
#