function yout=ie(f,Df,tspan,yin,options,varargin)
%IODE1 Solves ODEs using implicit nonadaptive method of order 1 (implicit Euler).
%  Y = IODE1(ODEFUN,JACOBIAN,TSPAN,Y0) integrates the system of ordinary
%  differential equations y' = f(t,y) using an implicit Euler method 
%  by stepping from T0 to T1 to T2, etc. Function ODEFUN(T,Y) must return 
%  f(t,y) in an N-dimensional column vector. The vector Y0 is the initial 
%  conditions at T0. Each row in the solution array Y corresponds to a time
%  specified in TSPAN. Function JACOBIAN(T,Y) must return the Jacobian 
%  matrix of f(t,y) which is needed to perform Newton iterations for 
%  solving a system of 2*N nonlinear algebraic equations each time step. 
%
%  Y = IODE1(ODEFUN,JACOBIAN,TSPAN,Y0,OPTIONS) solves as above with default
%  integration properties replaced by values in OPTIONS, structured array
%  where: 
%  
%  OPTIONS.tout specifies times at which the solution is to be returned.
%  This is meant either for stroboscopic sampling at integer multiples of 
%  TS or to return the solution at the final time specified in TSPAN. In 
%  the former case TSPAN = 0:TS/M:N*TS and tout = TS:TS:N*TS, whereas in
%  the latter tout = TSPAN(END). Default setting is tout = TSPAN.
%
%  OPTIONS.tol sets the error tolerance for the relative increment of 
%  the solution of the algebraic system using Newton's method. Default 
%  setting is tol = 1e-6.
%
%  OPTIONS.maxiter is the maximum number of iterations to be performed 
%  while solving the algebraic system. Default setting is maxiter = 15.
%
%  Y = IODE1(ODEFUN,JACOBIAN,TSPAN,Y0,OPTIONS,VARARGIN) passes the
%  additional parameters P1, P2, ... to ODEFUN and JACOBIAN as ODEFUN(T,Y,
%  P1,P2,...) and JACOBIAN(T,Y,P1,P2,...).
%
%  Example: Solve x'' + a^2 sin(x) = 0, x(0) = 2, x'(0) = 0, a = 1, using an
%           implicit Euler method with tolerance 1e-3 and 10 maximum Newton 
%           iterations per time step, and plot solution in phase space.
%
%           f=@(t,y,p) [y(2);-p^2*sin(y(1))];
%           Df=@(t,y,p) [0 1;-p^2*cos(y(1)) 0];
%           p=1;
%           y0=[2;0];
%           t=0:.025:4*pi;
%           options=struct('tout',t,'tol',1e-3,'maxiter',10);
%           y=iode1(f,Df,t,y0,options,p);
%           plot(y(:,1),y(:,2))

%  Francisco J. Beron-Vera, 2006/10/09

% Set default parameters.
if nargin<5
   tout=tspan;
   tol=1e-6;
   maxiter=15;
else
   tout=options.tout;
   tol=options.tol;
   maxiter=options.maxiter;
end

% Set dimensions.
Ny=length(yin);  
Ntout=length(tout);
Ntspan=length(tspan); 
if Ntspan==Ntout
   nsample=1;
   nmax=Ntout;
elseif (Ntspan>Ntout) && (Ntout>1)
   nsample=(Ntspan-1)/Ntout;
   nmax=nsample*Ntout+1;
else
   nsample=Ntspan-1;
   nmax=Ntspan;
end

% Fill-in solution array with zeros.
yout=zeros(Ntout,Ny); 
if Ntspan==Ntout
   yout(1,:)=yin'; 
   nout=2;
else
   nout=1;
end

% Main loop.
hspan=diff(tspan);
y=yin;
for n=2:nmax
   t=tspan(n-1);
   if ~mod(n,1000)
      disp(['t = ' num2str(t) ' of ' num2str(tspan(end))])
   end
   h=hspan(n-1);
   g=newton(@(g)  F(g, f,t,h,y,   varargin{:}), ...
            @(g) DF(g,Df,t,h,  Ny,varargin{:}), ...
            y,tol,maxiter);
   y=y+h*f(t+h,g,varargin{:});
   if ~mod(n-1,nsample)
      yout(nout,:)=y';
      nout=nout+1;
   end
end

%--------------------------------------------------------------------------
function out=F(g,f,t,h,y,varargin)
% F returns algebraic system F(g)=0.
out=g-y-h*f(t+h,g,varargin{:});

%--------------------------------------------------------------------------
function out=DF(g,Df,t,h,Ny,varargin)
% DF returns Jacobian matrix for F(g)=0.
out=eye(Ny)-h*Df(t+h,g,varargin{:});

%--------------------------------------------------------------------------
function xout=newton(sysfun,jacfun,xin,tol,maxiter,varargin)
% NEWTON solves nonlinear algebraic systems.
n=0;
dx=tol+1;
while (dx>tol) && (j<=maxiter)
   n=n+1;
   F=sysfun(xin,varargin{:});
   J=jacfun(xin,varargin{:});
   xout=xin-J\F;
   dx=norm(xout-xin);
   xin=xout;
end