# Implements the backwards Euler method using Newton iterations to solve for y_{n+1} in any the nonlinear function.
#
# An implementation of Newtons method:
#
#source('../../../../Conte/Code/Chapter3/utils.R')
source('http://www.waxworksmath.com/Authors/A_F/Conte/Code/Chapter3/utils.R')

beuler = function(ynp1_res_fn, ynp1_res_prime_fn, t0, y0, t1, h){
    #
    # Integrate from t0 -> t1 starting at y(t0) = y0 using a stepsize of h
    #
    n_max = ceiling(t1/h)
    ts = rep(NA, n_max)
    ys = rep(NA, n_max)
    
    tn = t0
    yn = y0
    for( n in 1:n_max ){
        # Use Newtons method to solve for y_{n+1}:
        #
        root_fn = function(y){ ynp1_res_fn(y, yn, tn, h) }
        root_fn_prime = function(y){ ynp1_res_prime_fn(y, yn, tn, h) }
        nm_res = newtons_method(yn, root_fn, root_fn_prime, xtol=1.e-3, ftol=1.e-3, ntol=10)
        
        ts[n] = t0 + n*h
        ys[n] = nm_res[1]
        
        # Prepare for the next call
        tn = ts[n]
        yn = ys[n]
    }

    # Put in the initial conditions into the array:
    ts = c( t0, ts )
    ys = c( y0, ys )
    list( t=ts, y=ys )
}