function phi = backward_euler( fname, phi0, t )
%phi = backward_euler( fname, phi0, t )
% 
% Integrates the ODE(s) whose RHS is given by fname forward in time.
%
% INPUTS:
%   fname - the name of the function to calculate the RHS.  It should take
%           two arguments: the time (t) and the solution value (phi).
%   phi0  - The initial solution value(s)
%   t     - The times at which the solution should be calculated.  The
%           timestep is determined as the difference between these times.
%
% OUTPUT:
%   phi - The solution value(s) at time tend.

% Author: James C. Sutherland
% Date:   Spring, 2007

nt       = length(t);  % number of time locations
[neq,nc] = size(phi0); % neq = number of equations.

if( nt<2 ) error('must specify at least two points for time.'); end
if( nc>1 ) error('initial guess must be a column vector'); end

% build an array to hold the solution at each time
phi = zeros(neq,nt);
phi(:,1) = phi0;

% set the error tolerance on the nonlinear solve.  Here we choose a
% tolerance that is relative to the initial condition of the solution.  If
% the initial solution is zero, then this reverts to 10^-9.
tol = 1e-9 * (1+phi0);

for i=2:nt

   % determine the time step
   dt = t(i)-t(i-1);
   
   % create the residual function that we will solve using Newton's method.
   % we create it inside this loop since it depends on dt, which is defined
   % inside the loop.
   f = @(y)( (y - phi0)/dt - feval(fname,t,y) );

   % solve the nonlinear system of equations for the solution variable at
   % the new time.
   phi(:,i) = newton_fd_sys( f, phi0, tol );
   
   % reset the guess for the next solution
   phi0 = phi(:,i);

end
