newtons_method = function(xn, root_fn, root_fn_prime, xtol=1.e-6, ftol=1.e-3, ntol=1000){
    # Algo 3.5: EPage 95
    #
    fn = root_fn(xn)
    xsize = abs(xn)
    fsize = abs(fn)
    
    n=0
    done = FALSE
    while( !done ){
        fn_prime = root_fn_prime(xn)
        stopifnot( abs(fn_prime)>1.e-9 ) # add an assert to flag if f'(x_n)=0
        xnp1 = xn - fn / fn_prime 

        # Check for convergence:
        #
        if( abs(xn-xnp1)<xtol*xsize || (abs(fn)/fsize)<ftol || n>=ntol ){
            done = TRUE
        }
        print(sprintf('Newtons method: n= %d; x_k= %f; x_kp1= %f', n, xn, xnp1))

        # Update root estimate:
        #
        xn = xnp1
        fn = root_fn(xn)

        n = n+1
    }

    return( c(xn, fn, n) )
}