function [tout,yout,varargout] = ode23t(odefile,tspan,y0,options,varargin)
%ODE23T  Solve moderately stiff ODEs and DAEs, trapezoidal rule.
%   [T,Y] = ODE23T('F',TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the
%   system of differential equations y' = F(t,y) from time T0 to TFINAL with
%   initial conditions Y0.  'F' is a string containing the name of an ODE
%   file.  Function F(T,Y) must return a column vector.  Each row in
%   solution array Y corresponds to a time returned in column vector T.  To
%   obtain solutions at specific times T0, T1, ..., TFINAL (all increasing
%   or all decreasing), use TSPAN = [T0 T1 ... TFINAL].
%   
%   [T,Y] = ODE23T('F',TSPAN,Y0,OPTIONS) solves as above with default
%   integration parameters replaced by values in OPTIONS, an argument
%   created with the ODESET function.  See ODESET for details.  Commonly
%   used options are scalar relative error tolerance 'RelTol' (1e-3 by
%   default) and vector of absolute error tolerances 'AbsTol' (all
%   components 1e-6 by default).
%   
%   [T,Y] = ODE23T('F',TSPAN,Y0,OPTIONS,P1,P2,...) passes the additional
%   parameters P1,P2,... to the ODE file as F(T,Y,FLAG,P1,P2,...) (see
%   ODEFILE).  Use OPTIONS = [] as a place holder if no options are set.
%   
%   It is possible to specify TSPAN, Y0 and OPTIONS in the ODE file (see
%   ODEFILE).  If TSPAN or Y0 is empty, then ODE23T calls the ODE file
%   [TSPAN,Y0,OPTIONS] = F([],[],'init') to obtain any values not supplied
%   in the ODE23T argument list.  Empty arguments at the end of the call
%   list may be omitted, e.g. ODE23T('F').
%   
%   The Jacobian matrix dF/dy is critical to reliability and efficiency.
%   Use ODESET to set JConstant 'on' if dF/dy is constant.  Set Vectorized
%   'on' if the ODE file is coded so that F(T,[Y1 Y2 ...]) returns
%   [F(T,Y1) F(T,Y2) ...].  Set JPattern 'on' if dF/dy is a sparse matrix
%   and the ODE file is coded so that F([],[],'jpattern') returns a sparsity
%   pattern matrix of 1's and 0's showing the nonzeros of dF/dy.  Set
%   Jacobian 'on' if the ODE file is coded so that F(T,Y,'jacobian') returns
%   dF/dy.
%   
%   As an example, the command
%   
%       ode23t('vdpode',[0 3000],[2 0],[],1000);
%   
%   solves the system y' = vdpode(t,y) with mu = 1000, using the default
%   relative error tolerance 1e-3 and the default absolute tolerance of 1e-6
%   for each component.  When called with no output arguments, as in this
%   example, ODE23T calls the default output function ODEPLOT to plot the
%   solution as it is computed.
%   
%   ODE23T can solve problems M(t,y)*y' = F(t,y) with a mass matrix M that
%   is nonsingular.  Use ODESET to set Mass to 'M', 'M(t)', or 'M(t,y)' if
%   the ODE file is coded so that F(T,Y,'mass') returns a constant,
%   time-dependent, or time- and state-dependent mass matrix, respectively.
%   The default value of Mass is 'none'.  See FEM1ODE, FEM2ODE, or
%   BATONODE.
%   
%   If M is singular, M(t)*y' = F(t,y) is a differential-algebraic equation
%   (DAE).  DAEs have solutions only when y0 is consistent, i.e. there is a
%   vector yp0 such that M(t0)*yp0 = f(t0,y0).  ODE23T can solve DAEs of
%   index 1 provided that M is not state dependent and y0 is sufficiently
%   close to being consistent.  You can use ODESET to set MassSingular to
%   'yes', 'no', or 'maybe'.  The default of 'maybe' causes ODE23T to test
%   whether the problem is a DAE. If it is, ODE23T treats y0 as a guess,
%   attempts to compute consistent initial conditions that are close to y0,
%   and goes on to solve the problem. When solving DAEs, it is very
%   advantageous to formulate the problem so that M is diagonal (a
%   semi-explicit DAE).  See HB1DAE or AMP1DAE.
%   
%   [T,Y,TE,YE,IE] = ODE23T('F',TSPAN,Y0,OPTIONS) with the Events property
%   in OPTIONS set to 'on', solves as above while also locating zero
%   crossings of an event function defined in the ODE file.  The ODE file
%   must be coded so that F(T,Y,'events') returns appropriate information.
%   See ODEFILE for details.  Output TE is a column vector of times at which
%   events occur, rows of YE are the corresponding solutions, and indices in
%   vector IE specify which event occurred.
%   
%   See also ODEFILE and
%       other ODE solvers:   ODE15S, ODE23S, ODE23TB, ODE45, ODE23, ODE113
%       options handling:    ODESET, ODEGET
%       output functions:    ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT
%       odefile examples:    VDPODE, BRUSSODE, B5ODE, CHM6ODE, FEM1ODE
%       Jacobian functions:  NUMJAC, COLGROUP

%   ODE23T is an implementation of the trapezoidal rule.  A "free"
%   interpolant is used. 

%   Mark W. Reichelt, Lawrence F. Shampine, Yanyuan Ma, and Jacek Kierzenka
%   12-18-97
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.9 $  $Date: 1998/05/18 15:06:37 $

true = 1;
false = ~true;

nsteps = 0;                             % stats
nfailed = 0;                            % stats
nfevals = 0;                            % stats
npds = 0;                               % stats
ndecomps = 0;                           % stats
nsolves = 0;                            % stats

if nargin == 0
  error('Not enough input arguments.  See ODE23T.');
elseif ~isstr(odefile) & ~isa(odefile, 'inline')
  error('First argument must be a single-quoted string.  See ODE23T.');
end

if nargin == 1
  tspan = []; y0 = []; options = [];
elseif nargin == 2
  y0 = []; options = [];
elseif nargin == 3
  options = [];
elseif ~isempty(options) & ~isa(options,'struct')
  if (length(tspan) == 1) & (length(y0) == 1) & (min(size(options)) == 1)
    tspan = [tspan; y0];
    y0 = options;
    options = [];
    varargin = {};
    msg = sprintf('Use ode23t(''%s'',tspan,y0,...) instead.',odefile);
    warning(['Obsolete syntax.  ' msg]);
  else
    error('Correct syntax is ode23t(''odefile'',tspan,y0,options).');
  end
end

% Get default tspan and y0 from odefile if none are specified.
if isempty(tspan) | isempty(y0)
  if (nargout(odefile) < 3) & (nargout(odefile) ~= -1)
    msg = sprintf('Use ode23t(''%s'',tspan,y0,...) instead.',odefile);
    error(['No default parameters in ' upper(odefile) '.  ' msg]);
  end
  [def_tspan,def_y0,def_options] = feval(odefile,[],[],'init',varargin{:});
  if isempty(tspan)
    tspan = def_tspan;
  end
  if isempty(y0)
    y0 = def_y0;
  end
  if isempty(options)
    options = def_options;
  else
    options = odeset(def_options,options);
  end
end

% Test that tspan is internally consistent.
tspan = tspan(:);
ntspan = length(tspan);
if ntspan == 1
  t0 = 0;
  next = 1;
else
  t0 = tspan(1);
  next = 2;
end
tfinal = tspan(ntspan);
if t0 == tfinal
  error('The last entry in tspan must be different from the first entry.');
end
tdir = sign(tfinal - t0);
if any(tdir * (tspan(2:ntspan) - tspan(1:ntspan-1)) <= 0)
  error('The entries in tspan must strictly increase or decrease.');
end

t = t0;
y = y0(:);
neq = length(y);
one2neq = (1:neq)';

% Get options, and set defaults.
rtol = odeget(options,'RelTol',1e-3);
if (length(rtol) ~= 1) | (rtol <= 0)
  error('RelTol must be a positive scalar.');
end
if rtol < 100 * eps 
  rtol = 100 * eps;
  warning(['RelTol has been increased to ' num2str(rtol) '.']);
end

atol = odeget(options,'AbsTol',1e-6);
if any(atol <= 0)
  error('AbsTol must be positive.');
end

normcontrol = strcmp(odeget(options,'NormControl','off'),'on');
if normcontrol
  if length(atol) ~= 1
    error('Solving with NormControl ''on'' requires a scalar AbsTol.');
  end
  normy = norm(y);
else
  if (length(atol) ~= 1) & (length(atol) ~= neq)
    error(sprintf(['Solving %s requires a scalar AbsTol, ' ...
                   'or a vector AbsTol of length %d'],upper(odefile),neq));
  end
  atol = atol(:);
end
threshold = atol / rtol;

% By default, hmax is 1/10 of the interval.
hmax = min(abs(tfinal-t), abs(odeget(options,'MaxStep',0.1*(tfinal-t))));
if hmax <= 0
  error('Option ''MaxStep'' must be greater than zero.');
end
htry = abs(odeget(options,'InitialStep'));
if htry <= 0
  error('Option ''InitialStep'' must be greater than zero.');
end

haveeventfun = strcmp(odeget(options,'Events','off'),'on');
if haveeventfun
  valt = feval(odefile,t,y,'events',varargin{:});
  teout = [];
  yeout = [];
  ieout = [];
end

if nargout > 0
  outfun = odeget(options,'OutputFcn');
else
  outfun = odeget(options,'OutputFcn','odeplot');
end
if isempty(outfun)
  haveoutfun = false;
else
  haveoutfun = true;
  outputs = odeget(options,'OutputSel',one2neq);
end
refine = odeget(options,'Refine',1);
printstats = strcmp(odeget(options,'Stats','off'),'on');

Janalytic = strcmp(odeget(options,'Jacobian','off'),'on');
Jconstant = strcmp(odeget(options,'JConstant','off'),'on');
vectorized = strcmp(odeget(options,'Vectorized','off'),'on');
Jpattern = strcmp(odeget(options,'JPattern','off'),'on');
if Jpattern
  Js = feval(odefile,[],[],'jpattern',varargin{:});
else
  Js = [];
end

mass = lower(odeget(options,'Mass','none'));

% Grandfathering code -- should be removed in a subsequent release.
Mconstant = odeget(options,'MassConstant');
if strcmp(mass,'on') | strcmp(mass,'off') | ...
      strcmp(Mconstant,'on') | strcmp(Mconstant,'off')
  if strcmp(mass,'on') | strcmp(mass,'off')
    warning(['Mass property values are ''none'', ''M'', ''M(t)'', ' ...
          'and ''M(t,y)'' (see ODESET).  Support for the''on'' and ' ...
          '''off'' values has been grandfathered and will disappear ' ...
          'in a future release.']);
  end
  if strcmp(Mconstant,'on') | strcmp(Mconstant,'off')
    warning(['The MassConstant property has been grandfathered and will ' ...
          'disappear in a future release (see ODESET).']);
  end
  if strcmp(Mconstant,'on')
    mass = 'm';
    warning('Assuming Mass value ''M'', a constant mass matrix.');
  elseif strcmp(mass,'on')
    mass = 'm(t)';
    warning('Assuming Mass value ''M(t)'', a time-dependent mass matrix.');
  else
    mass = 'none';
    warning('Assuming Mass value ''none'', no mass matrix.');
  end
end

switch(mass)
  case 'none', Mtype = 0;
  case 'm', Mtype = 1;
  case 'm(t)', Mtype = 2;
  case 'm(t,y)', Mtype = 3;
  otherwise, error('Unrecognized Mass property value.  See ODESET.');
end

Msingular = odeget(options,'MassSingular');
if Mtype == 3
  if strcmp(Msingular,'maybe')
    warning(['Solver assumes MassSingular is ''no'' for state-dependent ' ...
          'mass matrices.']);
    Msingular = 'no';
  elseif strcmp(Msingular,'yes')
    error(['MassSingular cannot be ''yes'' for state-dependent mass '...
          'matrices M(t,y).']);
  else
    Msingular = 'no';
  end
elseif isempty(Msingular)
  Msingular = 'maybe';
end

DAE = false;
RowScale = [];
SE = false;
if Mtype > 0
  Mt = feval(odefile,t,y,'mass',varargin{:});
  
  if nnz(Mt) == 0
    error('The mass matrix must have some non-zero entries.')
  end
  
  if ~strcmp(Msingular,'no')
    
    % If Mt is singular, DAE = true.  It is very advantageous to
    % to recognize a problem that is semi-explicit, SE = true.
    [r,c] = find(Mt);
    SE = isequal(r,c);                  % Test for a diagonal mass matrix.

    DAE = true;
    if strcmp(Msingular,'maybe')
      if SE
        DAE = any(diag(Mt) == 0);
      else
        DAE = (eps*nnz(Mt)*condest(Mt) > 1);
      end
    end
  end
  
else
  Mt = sparse(one2neq,one2neq,1,neq,neq); % Mtype = 0, mass matrix is sparse I.
end
Mcurrent = true;
Mtnew = Mt;

% Set the output flag.
if ntspan > 2
  outflag = 1;                          % output only at tspan points
elseif refine <= 1
  outflag = 2;                          % computed points, no refinement
else
  outflag = 3;                          % computed points, with refinement
  S = (1:refine-1)' / refine;
end

% Initialize method parameters.
% The first step is taken at order one with the backward Euler method (BDF1).
pow = 1/2;
gamma = 1;
start = true;

% The input arguments of odefile determine the args to use to evaluate f.
if (exist(odefile) == 3) | (nargin(odefile) == 2) % don't nargin MEX files
  args = {};                            % odefile accepts only (t,y)
else
  args = [{''} varargin];               % use (t,y,'',p1,p2,...)
end

f0 = feval(odefile,t,y,args{:});
nfevals = nfevals + 1;                  % stats
[m,n] = size(f0);
if n > 1
  error([upper(odefile) ' must return a column vector.'])
elseif m ~= neq
  msg = sprintf('an initial condition vector of length %d.',m);
  error(['Solving ' upper(odefile) ' requires ' msg]);
end

warnstat = warning;

% Compute the initial slope yp.  For DAEs the y input is usually just a
% guess.  Compute consistent initial conditions when necessary.
jthresh = atol + zeros(neq,1);
if DAE
   if SE | ~issparse(Mt)
     [y,yp,f0,dfdy,nFE,nPD,fac,g] = icsedae(odefile,t,SE,Mt,y,f0,rtol,...
         Janalytic,jthresh,vectorized,Js,varargin{:});
   else
     [y,yp,f0,dfdy,nFE,nPD,fac,g] = icdae(odefile,tspan,htry,Mt,y,f0,rtol,...
         Janalytic,jthresh,vectorized,Js,varargin{:});
  end
  nfevals = nfevals + nFE;
  npds = npds + nPD;
else
  if Mtype > 0
    [L,U] = lu(Mt);
    yp = U \ (L \ f0);
    ndecomps = ndecomps + 1;              % stats
    nsolves = nsolves + 1;                % stats
  else
    yp = f0;
  end
  
  if Janalytic
    dfdy = feval(odefile,t,y,'jacobian',varargin{:});
  else
    [dfdy,fac,g,nF] = ...
        numjac(odefile,t,y,f0,jthresh,[],vectorized,Js,[],args{:});
    nfevals = nfevals + nF;               % stats
  end
  npds = npds + 1;                        % stats
end
Jcurrent = true;
needNewJ = false;

% hmin is a small number such that t + hmin is clearly different from t in
% the working precision, but with this definition, it is 0 if t = 0.
hmin = 16*eps*abs(t);

if isempty(htry)
  % Compute an initial step size h using yp = y'(t).
  if normcontrol
    wt = max(normy,threshold);
    rh = 1.43 * (norm(yp) / wt) / rtol^pow;  % 1.43 = 1 / 0.7
  else  
    wt = max(abs(y),threshold);
    rh = 1.43 * norm(yp ./ wt,inf) / rtol^pow;
  end
  absh = min(hmax, abs(tspan(next) - t));
  if absh * rh > 1
    absh = 1 / rh;
  end
  absh = max(absh, hmin);
  
  if ~DAE
    % The error of BDF1 is 0.5*h^2*y''(t), so we can determine the optimal h.
    h = tdir * absh;
    tdel = (t + tdir*min(sqrt(eps)*max(abs(t),abs(t+h)),absh)) - t;
    f1 = feval(odefile,t+tdel,y,args{:});
    nfevals = nfevals + 1;              % stats
    dfdt = (f1 - f0) ./ tdel;
    if normcontrol
      if Mtype > 0
        rh = 1.43 * sqrt(0.5 * (norm(U \ (L \ (dfdt + dfdy*yp))) / wt) / rtol);
      else
        rh = 1.43 * sqrt(0.5 * (norm(dfdt + dfdy*yp) / wt) / rtol);
      end
    else
      if Mtype > 0
        rh = 1.43*sqrt(0.5*norm((U \ (L \ (dfdt+dfdy*yp))) ./ wt,inf) / rtol);
      else
        rh = 1.43 * sqrt(0.5 * norm((dfdt + dfdy*yp) ./ wt,inf) / rtol);
      end
    end
    absh = min(hmax, abs(tspan(next) - t));
    if absh * rh > 1
      absh = 1 / rh;
    end
    absh = max(absh, hmin);
  end
else
  absh = min(hmax, max(hmin, htry));
end
h = tdir * absh;

% Initialize the output function.
if haveoutfun
  feval(outfun,[t tfinal],y(outputs),'init');
end

% Allocate memory if we're generating output.
if nargout > 0
  if ntspan > 2                         % output only at tspan points
    tout = zeros(ntspan,1);
    yout = zeros(ntspan,neq);
  else                                  % alloc in chunks
    chunk = max(ceil(128 / neq),refine);
    tout = zeros(chunk,1);
    yout = zeros(chunk,neq);
  end
  nout = 1;
  tout(nout) = t;
  yout(nout,:) = y.';
end

% THE MAIN LOOP

z = h * yp;                             % z is the scaled derivative.
needNewLU = true;                       % Initialize LU.
done = false;
while ~done
  
  hmin = 16*eps*abs(t);
  abshlast = absh;
  absh = min(hmax, max(hmin, absh));
  h = tdir * absh;
  
  % Stretch the step if within 10% of tfinal-t.
  if 1.1*absh >= abs(tfinal - t)
    h = tfinal - t;
    absh = abs(h);
    done = true;
  end
  
  if absh ~= abshlast
    z = (absh / abshlast) * z;
    needNewLU = true;
  end
  
  % LOOP FOR ADVANCING ONE STEP.
  nofailed = true;                      % no failed attempts
  itfailed = 0;
  while true                            % Evaluate the formula.
    
    if normcontrol
      wt = max(normy,threshold);
    else
      wt = max(abs(y),threshold);
    end    
    
    if needNewJ
      if Janalytic
        dfdy = feval(odefile,t,y,'jacobian',varargin{:});
      else
        f0 = feval(odefile,t,y,args{:});
        [dfdy,fac,g,nF] = ...
            numjac(odefile,t,y,f0,jthresh,fac,vectorized,Js,g,args{:});
        nfevals = nfevals + nF + 1;     % stats
      end
      npds = npds + 1;                  % stats
      Jcurrent = true;
      needNewJ = false;
      if ~Mcurrent                      % possible only if state-dependent
        Mt = feval(odefile,t,y,'mass',varargin{:});
        Mcurrent = true;
      end
      Miter = Mt - (gamma*h)*dfdy;
      % Use explicit scaling of the equations when solving DAEs.
      if DAE
        RowScale = 1 ./ max(abs(Miter),[],2);
        Miter = sparse(one2neq,one2neq,RowScale) * Miter;
      end
      [L,U] = lu(Miter);
      ndecomps = ndecomps + 1;          % stats
      rate = [];
      needNewLU = false;
    elseif needNewLU
      if ~Mcurrent
        Mt = feval(odefile,t,y,'mass',varargin{:});
        Mcurrent = true;
      end
      Miter = Mt - (gamma*h)*dfdy;
      % Use explicit scaling of the equations when solving DAEs.
      if DAE
        RowScale = 1 ./ max(abs(Miter),[],2);
        Miter = sparse(one2neq,one2neq,RowScale) * Miter;
      end
      [L,U] = lu(Miter);
      ndecomps = ndecomps + 1;          % stats
      rate = [];
      needNewLU = false;
    end
    
    % Predict ynew, znew at tnew.
    tnew = t + h;
    if start                            % Backward Euler (BDF1)
      % Use linear interpolating polynomial.
      znew = z;
      ynew = y + z;
    else                                % Trapezoidal rule (TR)
      % Use quadratic interpolating polynomial.
      a1 = ((tnew - t_2)*(tnew - t_1))/((t - t_2)*(t - t_1));
      a2 = ((tnew - t_2)*(tnew - t))/((t_1 - t_2)*(t_1 - t));
      a3 = ((tnew - t_1)*(tnew - t))/((t_2 - t_1)*(t_2 - t));
      ynew = a1*y + a2*y_1 + a3*y_2;
      znew = 2*(ynew - y) - z;
    end
    
    % Mtnew is required in the RHS function evaluation.
    if Mtype == 2
      Mtnew = feval(odefile,tnew,ynew,'mass',varargin{:});
    end
    
    [ynew,znew,iter,itfail,rate] = ...
        itsolve(Mtnew,tnew,ynew,znew,gamma,h,L,U,odefile,rtol,wt,rate,...
                Mtype,DAE,RowScale,warnstat,varargin,args{:});
    nfevals = nfevals + iter;           % stats
    nsolves = nsolves + iter;           % stats
    
    if itfail
      nofailed = false;
      nfailed = nfailed + 1;            % stats
      itfailed = itfailed + 1;
      if Jcurrent                       % never false if Jconstant
        if absh <= hmin
          msg = sprintf(['Failure at t=%e.  Unable to meet integration ' ...
                         'tolerances without reducing the step size below ' ...
                         'the smallest value allowed (%e) at time t.\n'], ...
                        t,hmin);
          warning(msg);
          if haveoutfun
            feval(outfun,[],[],'done');
          end
          if printstats                 % print cost statistics
            fprintf('%g successful steps\n', nsteps);
            fprintf('%g failed attempts\n', nfailed);
            fprintf('%g function evaluations\n', nfevals);
            fprintf('%g partial derivatives\n', npds);
            fprintf('%g LU decompositions\n', ndecomps);
            fprintf('%g solutions of linear systems\n', nsolves);
          end
          if nargout > 0
            tout = tout(1:nout);
            yout = yout(1:nout,:);
            if haveeventfun
              varargout{1} = teout;
              varargout{2} = yeout;
              varargout{3} = ieout;
              varargout{4} = [nsteps;nfailed;nfevals; npds; ndecomps; nsolves];
            else
              varargout{1} = [nsteps;nfailed;nfevals; npds; ndecomps; nsolves];
            end
          end
          return;
          
        elseif itfailed == 3 & ~start
          start = true;                 % Need to restart.
          pow = 1/2;                    % Revert to BDF1.
          gamma = 1;
          
          % The local truncation error of BDF1 is (1/2)*y''*h^2. 
          % Approximate by differentiating a quadratic interpolant.
          d1 = y/((t-t_1)*(t-t_2));
          d2 = y_1/((t_1-t)*(t_1-t_2));
          d3 = y_2/((t_2-t)*(t_2-t_1));
          est = d1 + d2 + d3;
          err = norm(est ./ wt,inf)*absh^2;
          q = sqrt(err / rtol);
          ratio = hmax / absh;                 
          if 0.7 < q*ratio
            ratio = 0.7/q;
          end
          abshlast = absh;
          absh = max(ratio * absh, hmin);
          h = tdir * absh;
          z = (absh / abshlast) * z;    % Rescale z because of new h.
          needNewLU = true;
          done = false;  
        else
          abshlast = absh;
          absh = max(0.3 * absh, hmin);
          h = tdir * absh;
          z = (absh / abshlast) * z;    % Rescale z because of new h.
          needNewLU = true;
          done = false;
        end
      else   
        needNewJ = true;
      end
    else
      % Estimate the local truncation error.
      if normcontrol
        normynew = norm(ynew);
        wt = max(wt, normynew);
      else
        wt = max(wt, abs(ynew));
      end
      
      if start
        % The local truncation error is (1/2)*y''*h^2. 
        % Approximate by differencing the scaled derivative.
        err = norm(0.5*(znew - z) ./ wt,inf);
      else
        % The local truncation error is -(1/12)*y'''*h^3. 
        % Approximate by differentiating a cubic interpolant.
        c1 = ynew/((tnew-t_2)*(tnew-t_1)*(tnew-t));
        c2 = y/((t-tnew)*(t-t_2)*(t-t_1));
        c3 = y_1/((t_1-tnew)*(t_1-t)*(t_1-t_2));
        c4 = y_2/((t_2-t)*(t_2-t_1)*(t_2-tnew));
        est = (c1 + c2 + c3 + c4)/2;
        err = norm(est ./ wt,inf)*absh^3;
      end
      
      if err > rtol                     % Failed step
        nfailed = nfailed + 1;          % stats
        if absh <= hmin
          msg = sprintf(['Failure at t=%e.  Unable to meet integration ' ...
                         'tolerances without reducing the step size below ' ...
                         'the smallest value allowed (%e) at time t.\n'], ...
                        t,hmin);
          warning(msg);
          if haveoutfun
            feval(outfun,[],[],'done');
          end
          if printstats                 % print cost statistics
            fprintf('%g successful steps\n', nsteps);
            fprintf('%g failed attempts\n', nfailed);
            fprintf('%g function evaluations\n', nfevals);
            fprintf('%g partial derivatives\n', npds);
            fprintf('%g LU decompositions\n', ndecomps);
            fprintf('%g solutions of linear systems\n', nsolves);
          end
         if nargout > 0
            tout = tout(1:nout);
            yout = yout(1:nout,:);
            if haveeventfun
              varargout{1} = teout;
              varargout{2} = yeout;
              varargout{3} = ieout;
              varargout{4} = [nsteps;nfailed;nfevals; npds; ndecomps; nsolves];
            else
              varargout{1} = [nsteps;nfailed;nfevals; npds; ndecomps; nsolves];
            end
          end
          return;
        end
      
        nofailed = false;
        abshlast = absh;
        absh = max(abshlast * max(0.1, 0.7*(rtol/err)^pow), hmin);
        h = tdir * absh;
        z = (absh / abshlast) * z;
        needNewLU = true;
        done = false;
      else                              % Successful step
        break;
        
      end
    end

  end % while true
  nsteps = nsteps + 1;                  % stats
  
  if haveeventfun
    [te,ye,ie,valt,stop] = ...
        odezero('ntrp23t',odefile,valt,t,y,tnew,ynew,t0,varargin,h,z,znew);
    nte = length(te);
    if nte > 0
      if nargout > 2
        teout = [teout; te];
        yeout = [yeout; ye.'];
        ieout = [ieout; ie];
      end
      if stop                           % stop on a terminal event
        tnew = te(nte);
        ynew = ye(:,nte);
        done = true;
      end
    end
  end
  
  if nargout > 0
    oldnout = nout;
    if outflag == 2                     % computed points, no refinement
      nout = nout + 1;
      if nout > length(tout)
        tout = [tout; zeros(chunk,1)];
        yout = [yout; zeros(chunk,neq)];
      end
      tout(nout) = tnew;
      yout(nout,:) = ynew.';
    elseif outflag == 3                 % computed points, with refinement
      nout = nout + refine;
      if nout > length(tout)
        tout = [tout; zeros(chunk,1)];  % requires chunk >= refine
        yout = [yout; zeros(chunk,neq)];
      end
      i = oldnout+1:nout-1;
      tout(i) = t + (tnew-t)*S;
      yout(i,:) = ntrp23t(tout(i),t,y,tnew,ynew,h,z,znew).';
      tout(nout) = tnew;
      yout(nout,:) = ynew.';
    elseif outflag == 1                 % output only at tspan points
      while next <= ntspan
        if tdir * (tnew - tspan(next)) < 0
          if haveeventfun & done
            nout = nout + 1;
            tout(nout) = tnew;
            yout(nout,:) = ynew.';
          end
          break;
        elseif tnew == tspan(next)
          nout = nout + 1;
          tout(nout) = tnew;
          yout(nout,:) = ynew.';
          next = next + 1;
          break;
        end
        nout = nout + 1;                % tout and yout are already allocated
        tout(nout) = tspan(next);
        yout(nout,:) = ntrp23t(tspan(next),t,y,tnew,ynew,h,z,znew).';
        next = next + 1;
      end
    end
    
    if haveoutfun
      i = oldnout+1:nout;
      if ~isempty(i) & (feval(outfun,tout(i),yout(i,outputs).') == 1)
        tout = tout(1:nout);
        yout = yout(1:nout,:);
        if haveeventfun
          varargout{1} = teout;
          varargout{2} = yeout;
          varargout{3} = ieout;
          varargout{4} = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
        else
          varargout{1} = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
        end
        return;
      end
    end
    
  elseif haveoutfun
    if outflag == 2
      if feval(outfun,tnew,ynew(outputs)) == 1
        return;
      end
    elseif outflag == 3                 % computed points, with refinement
      tinterp = t + (tnew-t)*S;
      yinterp = ntrp23t(tinterp,t,y,tnew,ynew,h,z,znew);
      if feval(outfun,[tinterp; tnew],[yinterp(outputs,:), ynew(outputs)]) == 1
        return;
      end
    elseif outflag == 1                 % output only at tspan points
      ninterp = 0;
      while next <= ntspan 
        if tdir * (tnew - tspan(next)) < 0
          if haveeventfun & done
            ninterp = ninterp + 1;
            tinterp(ninterp,1) = tnew;
            yinterp(:,ninterp) = ynew;
          end
          break;
        elseif tnew == tspan(next)
          ninterp = ninterp + 1;
          tinterp(ninterp,1) = tnew;
          yinterp(:,ninterp) = ynew;
          next = next + 1;
          break;
        end
        ninterp = ninterp + 1;
        tinterp(ninterp,1) = tspan(next);
        yinterp(:,ninterp) = ntrp23t(tspan(next),t,y,tnew,ynew,h,z,znew);
        next = next + 1;
      end
      if ninterp > 0
        if feval(outfun,tinterp(1:ninterp),yinterp(outputs,1:ninterp)) == 1
          return;
        end
      end
    end
  end
  
  % Advance the integration one step.
  if start
    % Generate "previously computed" solution by interpolation and
    % put the stored solutions in order, though this is not necessary.
    t_1 = t + 0.5*(tnew - t);
    y_1 = ntrp23t(t_1,t,y,tnew,ynew,h,z,znew);
    t_2 = t;
    y_2 = y;
    
    t = tnew;
    y = ynew;
    if normcontrol
      normy = normynew;
    end
    z = znew;
    Jcurrent = Jconstant;
    switch Mtype
    case {0,1}
      Mcurrent = true;                  % Constant mass matrix I or M.
    case 2
      % M(t) has already been evaluated at tnew in Mtnew.
      Mt = Mtnew;
      Mcurrent = true;
    case 3
      % M(t,y) has not yet been evaluated at the accepted ynew.
      Mcurrent = false;
    end
    
    % Change from BDF1 to TR.  Doubling the step size is reasonable 
    % with a raise in order and does not change the LU decomposition. 
    start = false;
    pow = 1/3;
    gamma = 0.5;
    absh = 2 * absh;
    z = 2 * z;
  else
    t_2 = t_1;
    y_2 = y_1;
    t_1 = t;
    y_1 = y;
    
    t = tnew;
    y = ynew;
    if normcontrol
      normy = normynew;
    end
    z = znew;
    Jcurrent = Jconstant;
    switch Mtype
    case {0,1}
      Mcurrent = true;                  % Constant mass matrix I or M.
    case 2
      % M(t) has already been evaluated at tnew in Mtnew.
      Mt = Mtnew;
      Mcurrent = true;
    case 3
      % M(t,y) has not yet been evaluated at the accepted ynew.
      Mcurrent = false;
    end
    
    if nofailed
      q = (err / rtol)^pow;
      ratio = hmax / absh;                 
      if 0.7 < q*ratio
        ratio = 0.7/q;
      end
      ratio = min(5, max(0.2, ratio));
      if abs(ratio - 1) > 0.2
        absh = ratio * absh;
        needNewLU = true;
        z = ratio * z;
      end
    end
  end
  
end % while ~done

if haveoutfun
  feval(outfun,[],[],'done');
end

if printstats                           % print cost statistics
  fprintf('%g successful steps\n', nsteps);
  fprintf('%g failed attempts\n', nfailed);
  fprintf('%g function evaluations\n', nfevals);
  fprintf('%g partial derivatives\n', npds);
  fprintf('%g LU decompositions\n', ndecomps);
  fprintf('%g solutions of linear systems\n', nsolves);
end

if nargout > 0
  tout = tout(1:nout);
  yout = yout(1:nout,:);
  if haveeventfun
    varargout{1} = teout;
    varargout{2} = yeout;
    varargout{3} = ieout;
    varargout{4} = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
  else
    varargout{1} = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
  end
end

%------------------------------------------------------------------------------

function [y,z,iter,itfail,rate] = ...
    itsolve(M,t,y,z,gamma,h,L,U,odefile,rtol,wt,rate,...
            Mtype,DAE,RowScale,warnstat,args,varargin)

% Solve the nonlinear equation M*z = h*f(t,v+gamma*z) and y = v+gamma*z.  
% The value v is incorporated in the predicted y and is not needed because
% the y is corrected using corrections to z.  The argument t is constant 
% during the iteration.  The function f(t,y) is evaluated by odefile.m. 
% L,U is the lu decomposition of the matrix M - gamma*h*dfdy, where dfdy
% is an approximate Jacobian of f.  A simplified Newton (chord) iteration 
% is used, so dfdy and the decomposition are held constant.  z is computed 
% to an accuracy of kappa*rtol.  The rate of convergence of the iteration 
% is estimated. If the iteration succeeds, itfail is set false and the
% estimated rate is returned for use on a subsequent step.  rate can be
% used as long as neither h nor dfdy changes. 

maxiter = 5;
kappa = 0.5;
itfail = 0;
minnrm = 100 * eps * norm(y ./ wt,inf);

for iter = 1:maxiter
  if Mtype == 3
    M = feval(odefile,t,y,'mass',args{:});
  end
  rhs = h * feval(odefile,t,y,varargin{:}) - M * z;
  if DAE                                % Account for row scaling.
    rhs = RowScale .* rhs;
  end
  warning('off');
  del = U \ (L \ rhs);
  warning(warnstat);
  z = z + del;
  y = y + gamma*del;
  newnrm = norm(del ./ max(wt,abs(y)),inf);
  
  if newnrm <= minnrm
    break;
  elseif iter == 1
    if ~isempty(rate)
      errit = newnrm * rate / (1 - rate) ;
      if errit <= 0.1*kappa*rtol
        break;
      end
    else
      rate = 0;
    end
  elseif newnrm > 0.9*oldnrm
    itfail = 1;
    break;
  else
    rate = max(0.9*rate, newnrm / oldnrm);
    errit = newnrm * rate / (1 - rate);
    if errit <= kappa*rtol
      break;
    elseif iter == maxiter
      itfail = 1;
      break;
    elseif kappa*rtol < errit*rate^(maxiter-iter)
      itfail = 1;
      break;
    end
  end
  
  oldnrm = newnrm;
end
