function [t,x] = ieuler(Ffile,DFfile,t0,tf,x0,h)
% function [t,x] = ieuler(Ffile,DFfile,t0,tf,x0,h)               4/23/92
% Integrates ODE's backwards with the Implicit Euler Method
% (implicit equations are solved via linear approximation)
%
% For the system of ODE's:   dx/dt = F(x)
% IEULER.M iterates:         x[k+1] = x[k] - h [I + h DF(x[k])]^(-1) F(x[k])
% starting from x[0] = x0  (x is a column vector)
% See Mulligan, Casey B.  "Progressive Paths" equation (4-10)
%
% INPUTS:
% Ffile  = filename of F(x), a string
% DFfile = filename of DF(x), a string
% Ffile and DFfile should be function files with one (vector) argument, and
% Ffile should have column vector output of the same dimension as x
% DFfile should output a square matrix, a row the partials of one component
% Integration is from t0 to tf, at step size h
% OUTPUTS:
% t is a row vector of times [t0 (t0+h) (t0+2h) ... tf]
% x is a matrix of points on the trajectory [x(t0) x(t0+h) x(t0+2h) ...x(tf)]

xk = x0; tk = t0; x = xk; t = t0; n = length(x0);
while tk <= tf
     eval(['Fofx = ' Ffile '(xk);']);
     eval(['DFofx = ' DFfile '(xk);']);
     xk = xk - (eye(n) + h*DFofx)\Fofx*h;
     tk = tk + h;
     t = [t tk];  x = [x xk];
end
