#
# Borrowed and modified from: http://www.theresearchkitchen.com/archives/679
#
improved_euler = 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
        k1 = dy.dx(x, ys[i])
        k2 = dy.dx(x+h, ys[i] + h*k1)
        ys[i+1] = ys[i] + (h/2)*(k1 + k2)
        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 = improved_euler(dy.dx, start=0, end=0.5, h=0.1, y0=3)
#
# res$ys
# 
# [1] 3.000000 3.490000 3.961950 4.417565 4.858396 5.285848
#