% One step integration of an ODE by Euler's method
% CCK 2009-02-12
function x1 = backward_euler(ode, x0, t0, h, J)
% f = function handle for ODE,  dx/dt = f(x,t)
% x0 = initial x
% t0 = initial time
% h = timestep (t1 - t0)
% J = jacobian matrix (in the same format expected by Octave's LSODE).
%     If this argument is used, then J is put into a persistent variable
%     and the rest of the calculation is aborted.
%
% The result x1 is x(t1).


persistent Jacobian
if exist('J','var')
   Jacobian = J;
   x1 = 'The Jacobian has been set.';
   return
end

N = numel(x0); % Size of state vector
t1 = t0 + h;

% Implicit time step:
x1 = x0 + h * ( ( eye(N) - h*Jacobian(x0,t1) ) \ (ode(x0,t1)') )';
   % The \ operator inverts the Jacobian and performs matrix multiplication.

% Note that ode() had to be turned into a column vector, and then I had to
% convert the result of the multiplication back into a row vector to be
% compatible with the state vector. I begin to see why Octave's LSODE
% implementation is standardized around column vectors.
