source('runge_kutta.R')

predictor_corrector = 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)
    fs = numeric(nsteps+1)
    xs[1] = start
    ys[1] = y0
    fs[1] = dy.dx(xs[1], ys[1])

    # Get the first 4 values from the Runge-Kutta method
    #
    rk = runge_kutta(dy.dx, h, y0, start, end=start+4*h)
    for (i in 2:5){
        xs[i] = start + (i-1)*h
        ys[i] = rk$ys[i]
        fs[i] = dy.dx(xs[i], ys[i])
    }

    # Get the next values using a predictive corrector method:
    #
    for (i in 5:nsteps){
        x = start + (i-1)*h
        xs[i+1] = x + h

        # Take a 4th order Adams-Bashforth step:
        #
        y_predict = ys[i] + (h/24)*(55*fs[i] - 59*fs[i-1] + 37*fs[i-2] - 9*fs[i-3])

        # Correct this with a 4th order Adams-Moulton step: 
        # 
        f_predict = dy.dx(x+h, y_predict)
        ys[i+1] = ys[i] + (h/24)*(9*f_predict + 19*fs[i] - 5*fs[i-1] + fs[i-2])

        fs[i+1] = dy.dx(x + h, ys[i+1])
    }
    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 = predictor_corrector(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.286939
#