vector_predictor_corrector = function(vector_f=function(t,x,y){}, h=1E-7, xy0=c(0, 0), start=0, end=1){
    #
    # xy0 = initial condition [x(t_0), y(t_0)]
    # start = t_0 = start of time interation
    # end = t_1 = end of time integration
    # h = stepsize
    #
    nsteps = (end-start)/h
    ts = rep(NA, nsteps+1)
    xys = matrix(rep(NA, 2*(nsteps+1)), nrow=2, ncol=nsteps+1) # (2=n_states) x n_timesteps
    fs = matrix(rep(NA, 2*(nsteps+1)), nrow=2, ncol=nsteps+1) # (2=n_states) x n_timesteps
    ts[1] = start
    xys[1, 1] = xy0[1]
    xys[2, 1] = xy0[2]
    fs[, 1] = vector_f(ts[1], xys[1, 1], xys[2, 1])

    # Get the first 4 values from Euler's method (better methods could be used): 
    #
    for (i in 2:5){
        ts[i] = start + (i-1)*h
        f = vector_f(ts[i-1], xys[1, i-1], xys[2, i-1])
        xys[1, i] = xys[1, i-1] + h * f[1]
        xys[2, i] = xys[2, i-1] + h * f[2]
        fs[, i] = vector_f(ts[i], xys[1, i], xys[2, i])
    }

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

        f = vector_f(ts[i-1], xys[1, i-1], xys[2, i-1]) # this is f_n

        # Take a 4th order Adams-Bashforth step:
        #
        xy_predict = xys[, 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 = vector_f(ts[i+1], xy_predict[1], xy_predict[2])
        xys[, i+1] = xys[, i] + (h/24)*(9*f_predict + 19*fs[, i] - 5*fs[, i-1] + fs[, i-2])
        
        fs[, i+1] = vector_f(ts[i+1], xys[1, i+1], xys[2, i+1])
    } 
    list( ts=ts, xys=xys, fs=fs )
}