function [ tv, uv ] = backward_euler( fun, dfun_y, y0, t0, tf, Nh )
% BACKWARD_EULER Backward Euler method for the scalar ODE in the form
% y'(t) = f(t,y(t)),  t \in (t0,tf)
% y(0) = y_0
%
% The Newton method is used to solve the nonlinear equation at each time
% step. The function newton.m is used.
% 
%  [ tv, uv ] = backward_euler( fun, dfun_y, y0, t0, tf, Nh )
%  Inputs: fun    = function handle for f(t,y), fun = @(t,y) ...
%          dfun_y = derivative of f(t,y) w.r.t. y, dfun_y = @(t,y) ...
%          y0     = initial value
%          t0     = initial time
%          tf     = final time
%          Nh     = number of time subintervals
%  Output: tv     = vector of time steps (1 x (Nh+1))
%          uv     = vector of approximate solution at times tv
%

tv = linspace( t0, tf, Nh + 1 );
h = ( tf - t0 ) / Nh;

uv = zeros( 1, Nh + 1 );
uv( 1 ) = y0;

for n = 1 : Nh    
    f = @( x ) x - uv( n ) - h * fun( tv( n + 1 ), x );
    df = @( x ) 1 - h * dfun_y( tv( n + 1 ), x );
    [ xv, res, niter ] = newton( f, df, uv( n ), 1e-10, 20 );            
    uv( n + 1 ) = xv( end );
end

return