function varargout = ode23tb(ode,tspan,y0,options,varargin)
%ODE23TB Solve stiff differential equations, low order method.
%   [TOUT,YOUT] = ODE23TB(ODEFUN,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. ODEFUN is a function handle. For a scalar T 
%   and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
%   to f(t,y). Each row in the solution array YOUT corresponds to a time
%   returned in the column vector TOUT.  To obtain solutions at specific
%   times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
%   [T0 T1 ... TFINAL].      
%   
%   [TOUT,YOUT] = ODE23TB(ODEFUN,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).
%   If certain components of the solution must be non-negative, use
%   ODESET to set the 'NonNegative' property to the indices of these
%   components.  The 'NonNegative' property is ignored for problems
%   where there is a mass matrix. 
%      
%   The Jacobian matrix df/dy is critical to reliability and efficiency. Use
%   ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns 
%   the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. 
%   If the 'Jacobian'option is not set (the default), df/dy is approximated 
%   by finite differences. Set 'Vectorized' 'on' if the ODE function is coded 
%   so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. 
%   If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of
%   df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y)
%   depends on component j of y, and 0 otherwise.    
%   
%   ODE23TB can solve problems M(t,y)*y' = f(t,y) with mass matrix M(t,y)
%   that is nonsingular. Use ODESET to set the 'Mass' property to a function
%   handle MASS if MASS(T,Y) returns the value of the mass matrix. If the mass
%   matrix is constant, the matrix can be used as the value of the 'Mass'
%   option. Problems with state-dependent mass matrices are more
%   difficult. If the mass matrix does not depend on the state variable Y and
%   the function MASS is to be called with one input argument T, set
%   'MStateDependence' to 'none'. If the mass matrix depends weakly on Y, set
%   'MStateDependence' to 'weak' (the default) and otherwise, to 'strong'. In
%   either case the function MASS is to be called with the two arguments
%   (T,Y). If there are many differential equations, it is important to
%   exploit sparsity: Return a sparse M(t,y). Either supply the sparsity
%   pattern of df/dy using the 'JPattern' property or a sparse df/dy using
%   the Jacobian property. For strongly state-dependent M(t,y), set
%   'MvPattern' to a sparse matrix S with S(i,j) = 1 if for any k, the (i,k)
%   component of M(t,y) depends on component j of y, and 0 otherwise. ODE15S
%   and ODE23T can solve problems with singular mass matrices. 
%   
%   [TOUT,YOUT,TE,YE,IE] = ODE23TB(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
%   property in OPTIONS set to a function handle EVENTS, solves as above 
%   while also finding where functions of (T,Y), called event functions, 
%   are zero. For each function you specify whether the integration is 
%   to terminate at a zero and whether the direction of the zero crossing 
%   matters. These are the three column vectors returned by EVENTS: 
%   [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
%   VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
%   is to terminate at a zero of this event function and 0 otherwise. 
%   DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
%   zeros where the event function is increasing, and -1 if only zeros where 
%   the event function is decreasing. 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.        
%
%   SOL = ODE23TB(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
%   used with DEVAL to evaluate the solution or its first derivative at 
%   any point between T0 and TFINAL. The steps chosen by ODE23TB are returned 
%   in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
%   the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
%   of points at which events occurred. Columns of SOL.ye are the corresponding 
%   solutions, and indices in vector SOL.ie specify which event occurred. 
%   
%   Example
%         [t,y]=ode23tb(@vdp1000,[0 3000],[2 0]);   
%         plot(t,y(:,1));
%     solves the system y' = vdp1000(t,y), using the default relative error
%     tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
%     component, and plots the first component of the solution.
%
%   See also ODE15S, ODE23S, ODE23T, ODE45, ODE23, ODE113, ODE15I,
%            ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
%            ODEEXAMPLES, VDPODE, FEM1ODE, BRUSSODE, FUNCTION_HANDLE.

%   ODE23TB is an implementation of TR-BDF2, an implicit Runge-Kutta 
%   formula with a first stage that is a trapezoidal rule (TR) step and 
%   a second stage that is a backward differentiation formula (BDF) of 
%   order two.  By construction, the same iteration matrix is used in 
%   evaluating both stages.  The formula was proposed by Bank and Rose.
%   Here the improved error estimator and more efficient evaluation of
%   M. Hosea and L.F. Shampine are used.  A "free" interpolant is used.

%   Mark W. Reichelt, Lawrence F. Shampine, and Yanyuan Ma, 7-1-97
%   Copyright 1984-2023 The MathWorks, Inc.

solver_name = 'ode23tb';

if nargin < 4
    options = [];
    if nargin < 3
        y0 = [];
        if nargin < 2
            tspan = [];
        end
    end
end

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

[ode, odeIsFuncHandle, odeTreatAsMFile] = packageAsFuncHandle(ode);

% Output
output_sol = (odeIsFuncHandle && (nargout==1));      % sol = odeXX(...)
output_ty  = (~output_sol && (nargout > 0));  % [t,y,...] = odeXX(...)
% There might be no output requested...

sol = []; t2data = []; y2data = [];
if output_sol
    sol.solver = solver_name;
    sol.extdata.odefun = ode;
    sol.extdata.options = options;
    sol.extdata.varargin = varargin;
end

% Handle solver arguments
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...
    options, threshold, rtol, normcontrol, normy, hmax, htry, htspan] =  ...
    odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
nfevals = nfevals + 1;

% Handle the output
if nargout > 0
    outputFcn = odeget(options,'OutputFcn',[]);
else
    outputFcn = odeget(options,'OutputFcn',@odeplot);
end
outputArgs = {};
if isempty(outputFcn)
    haveOutputFcn = false;
else
    haveOutputFcn = true;
    outputs = odeget(options,'OutputSel',1:neq);
    if isa(outputFcn,'function_handle')
        % With MATLAB 6 syntax pass additional input arguments to outputFcn.
        outputArgs = varargin;
    end
end
refine = max(1,odeget(options,'Refine',1));
if ntspan > 2
    outputAt = 'RequestedPoints';         % output only at tspan points
elseif refine <= 1
    outputAt = 'SolverSteps';             % computed points, no refinement
else
    outputAt = 'RefinedSteps';            % computed points, with refinement
    S = (1:refine-1) / refine;
end
printstats = odeget(options,'Stats','off') == "on";

% Handle the event function
[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout] = ...
    odeevents(odeIsFuncHandle,ode,t0,y0,options,varargin);

% Handle the mass matrix
[Mtype, Mt, Mfun, Margs, dMoptions] = odemass(odeIsFuncHandle,ode,t0,y0,...
    options,varargin);
if Mtype > 0
    Msingular = odeget(options,'MassSingular','no');
    if strcmp(Msingular,'maybe')
        warning(message('MATLAB:ode23tb:MassSingularAssumedNo'));
    elseif strcmp(Msingular,'yes')
        error(message('MATLAB:ode23tb:MassSingularYes'));
    end
end

% Non-negative solution components
idxNonNegative = odeget(options,'NonNegative',[]);
nonNegative = ~isempty(idxNonNegative);
if nonNegative
    if Mtype == 0
        % Explicit ODE -- modify the derivative function
        [ode,thresholdNonNegative] = odenonnegative(ode,y0,threshold,idxNonNegative);
        f0 = ode(t0,y0,odeArgs{:});
        nfevals = nfevals + 1;
    else
        % Linearly implicit ODE/DAE -- ignore non-negativity constraints
        warning(message('MATLAB:ode23tb:NonNegativeIgnoredForLinearlyImplicitSystems'));
        nonNegative = false;
        idxNonNegative = [];
    end
end

% Handle the Jacobian
[Jconstant,Jac,Jargs,Joptions] = ...
    odejacobian(odeIsFuncHandle,ode,t0,y0,options,varargin);
Janalytic = isempty(Joptions);

% if not set via 'options', initialize constant Jacobian here
if Jconstant
    if isempty(Jac) % use odenumjac
        [Jac,Joptions.fac,nF] = odenumjac(ode, {t0,y0,odeArgs{:}}, f0, Joptions); %#ok<CCAT>
        nfevals = nfevals + nF;
        npds = npds + 1;
    elseif ~isa(Jac,'numeric')  % not been set via 'options'
        Jac = Jac(t0,y0,Jargs{:}); % replace by its value
        npds = npds + 1;
    end
end

t = t0;
y = y0;

Mcurrent = true;
needNewM = false;
Mt2 = Mt;
Mtnew = Mt;

% Initialize method parameters.
pow = 1/3;
alpha = 2 - sqrt(2);
d = alpha/2;
gg = sqrt(2)/4;

% Coefficients of the error estimate.
c1 = (alpha - 1)/3;
c2 = 1/3;
c3 = -alpha/3;

% Coefficients for the predictors
p31 = 1.5 + sqrt(2);
p32 = 2.5 + 2*sqrt(2);
p33 = - (6 + 4.5*sqrt(2));

% Compute the initial slope yp.
if Mtype > 0
    [Factors, piv] = odefactorize(Mt);
    yp = odesolve(Factors, piv, f0);
    ndecomps = ndecomps + 1;
    nsolves = nsolves + 1;
else
    yp = f0;
end

if Jconstant
    dfdy = Jac;
elseif Janalytic
    dfdy = Jac(t,y,Jargs{:});
    npds = npds + 1;
else
    [dfdy,Joptions.fac,nF] = odenumjac(ode, {t,y,odeArgs{:}}, f0, Joptions); %#ok<CCAT>
    nfevals = nfevals + nF;
    npds = npds + 1;
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, htspan);
    if absh * rh > 1
        absh = 1 / rh;
    end
    absh = max(absh, hmin);

    % Estimate error of first order Taylor series, 0.5*h^2*y''(t),
    % and use rule of thumb to select step size for second order method.
    h = tdir * absh;
    tdel = (t + tdir*min(sqrt(eps)*max(abs(t),abs(t+h)),absh)) - t;
    f1 = ode(t+tdel,y,odeArgs{:});
    nfevals = nfevals + 1;
    dfdt = (f1 - f0) ./ tdel;
    DfDt = dfdt + dfdy*yp;
    if normcontrol
        if Mtype > 0
            rh = 1.43*sqrt(0.5 * (norm(odesolve(Factors, piv, DfDt)) / wt)) / rtol^pow;
        else
            rh = 1.43 * sqrt(0.5 * (norm(DfDt) / wt)) / rtol^pow;
        end
    else
        if Mtype > 0
            rh = 1.43*sqrt(0.5*norm(odesolve(Factors, piv, DfDt)./ wt,inf)) / rtol^pow;
        else
            rh = 1.43 * sqrt(0.5 * norm( DfDt ./ wt,inf)) / rtol^pow;
        end
    end
    absh = min(hmax, htspan);
    if absh * rh > 1
        absh = 1 / rh;
    end
    absh = max(absh, hmin);
else
    absh = min(hmax, max(hmin, htry));
end
h = tdir * absh;

% Allocate memory if we're generating output.
nout = 0;
tout = []; yout = [];
if nargout > 0
    if output_sol
        chunk = min(max(100,50*refine), refine+floor((2^12)/neq));
        tout = zeros(1,chunk);
        yout = zeros(neq,chunk);
        t2data = zeros(1,chunk);
        y2data = zeros(neq,chunk);
    else
        if ntspan > 2                         % output only at tspan points
            tout = zeros(1,ntspan);
            yout = zeros(neq,ntspan);
        else                                  % alloc in chunks
            chunk = min(max(100,50*refine), refine+floor((2^13)/neq));
            tout = zeros(1,chunk);
            yout = zeros(neq,chunk);
        end
    end
    nout = 1;
    tout(nout) = t;
    yout(:,nout) = y;
end

% Initialize the output function.
if haveOutputFcn
    feval(outputFcn,[t tfinal],y(outputs),'init',outputArgs{:});
end

% THE MAIN LOOP

z = h * yp;                             % z is the scaled derivative.
if Mtype == 4
    [dMzdy,dMoptions.fac] = odenumjac(@odemxv, {Mfun,t,y,z,Margs{:}}, Mt*z, ...
        dMoptions); %#ok<CCAT>
end
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
    while true                            % Evaluate the formula.

        if normcontrol
            wt = max(normy,threshold);
        else
            wt = max(abs(y),threshold);
        end

        if needNewJ
            if Janalytic
                dfdy = Jac(t,y,Jargs{:});
            else
                f0 = ode(t,y,odeArgs{:});
                [dfdy,Joptions.fac,nF] = odenumjac(ode, {t,y,odeArgs{:}}, f0, Joptions); %#ok<CCAT>
                nfevals = nfevals + nF + 1;
            end
            npds = npds + 1;
            Jcurrent = true;
            needNewJ = false;
            needNewLU = true;
        end
        if needNewM
            Mt = Mfun(t,y,Margs{:});
            Mcurrent = true;
            if Mtype == 4
                [dMzdy,dMoptions.fac] = odenumjac(@odemxv, {Mfun,t,y,z,Margs{:}}, Mt*z, ...
                    dMoptions); %#ok<CCAT>
            end
            needNewM = false;
            needNewLU = true;
        end
        if needNewLU
            Miter = Mt - (d*h)*dfdy;
            if Mtype == 4
                Miter = Miter + dMzdy;
            end
            [Factors, piv] = odefactorize(Miter);
            ndecomps = ndecomps + 1;
            rate = [];
            needNewLU = false;
        end

        % The first stage is a TR step from t to t2.
        t2 = t + alpha*h;
        y2 = y + alpha*z;
        z2 = z;

        % Mt2 is required in the RHS function evaluation.
        if Mtype == 2   % M(t)
            if odeIsFuncHandle
                Mt2 = Mfun(t2,Margs{:});
            else
                Mt2 = Mfun(t2,y2,Margs{:});
            end
        end

        [y2,z2,iter,itfail1,rate] = ...
            itsolve(Mt2,t2,y2,z2,d,h,Factors,piv,ode,...
            rtol,wt,rate,Mtype,Mfun,Margs,odeArgs);
        nfevals = nfevals + iter;
        nsolves = nsolves + iter;
        itfail2 = false;                    % make sure well-defined later
        if ~itfail1
            % The second stage is a step from t2 to tnew with BDF2.
            if normcontrol
                wt = max(wt,norm(y2));
            else
                wt = max(wt,abs(y2));
            end

            tnew = t + h;
            if done
                tnew = tfinal;   % Hit end point exactly.
            end
            znew = p31*z + p32*z2 + p33*(y2 - y);
            ynew = y + gg * (z + z2) + d * znew;

            % Mtnew is required in the RHS function evaluation.
            if Mtype == 2   % M(t)
                if odeIsFuncHandle
                    Mtnew = Mfun(tnew,Margs{:});
                else
                    Mtnew = Mfun(tnew,ynew,Margs{:});
                end
            end

            [ynew,znew,iter,itfail2,rate] = ...
                itsolve(Mtnew,tnew,ynew,znew,d,h,Factors,piv,ode,rtol, ...
                wt,rate,Mtype,Mfun,Margs,odeArgs);
            nfevals = nfevals + iter;
            nsolves = nsolves + iter;
        end

        if itfail1 || itfail2                % Unable to evaluate a stage.
            nofailed = false;
            nfailed = nfailed + 1;
            if Jcurrent && Mcurrent
                if absh <= hmin
                    warning(message('MATLAB:ode23tb:IntegrationTolNotMet', sprintf( '%e', t ), sprintf( '%e', hmin )));
                    solver_output = odefinalize(solver_name, sol,...
                        outputFcn, outputArgs,...
                        printstats, [nsteps, nfailed, nfevals,...
                        npds, ndecomps, nsolves],...
                        nout, tout, yout,...
                        haveEventFcn, teout, yeout, ieout,...
                        {t2data,y2data,idxNonNegative},...
                        t);
                    if nargout > 0
                        varargout = solver_output;
                    end
                    return;
                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 = ~Jcurrent;
                needNewM = ~Mcurrent;
            end
        else
            % Estimate the local truncation error.
            if normcontrol
                normynew = norm(ynew);
                wt = max(wt, normynew);
            else
                wt = max(wt, abs(ynew));
            end

            est1 = c1*z + c2*z2 + c3*znew;
            err1 = norm(est1 ./ wt,inf);
            % Modify the estimate to improve it at infinity.  With this a
            % larger step size, but not "too" much larger, is reasonable.
            est2 = odesolve(Factors, piv, est1);
            nsolves = nsolves + 1;
            err2 = norm(est2 ./ wt,inf);
            err = max(err2, err1 / 16);

            NNrejectStep = false;
            if nonNegative && (err <= rtol) && any(ynew(idxNonNegative)<0)
                if normcontrol
                    errNN = norm( max(0,-ynew(idxNonNegative)) / wt );
                else
                    errNN = norm( max(0,-ynew(idxNonNegative)) ./ thresholdNonNegative, inf);
                end
                if errNN > rtol
                    err = errNN;
                    NNrejectStep = true;
                end
            end

            if err > rtol                     % Failed step
                nfailed = nfailed + 1;
                if absh <= hmin
                    warning(message('MATLAB:ode23tb:IntegrationTolNotMet', sprintf( '%e', t ), sprintf( '%e', hmin )));
                    solver_output = odefinalize(solver_name, sol,...
                        outputFcn, outputArgs,...
                        printstats, [nsteps, nfailed, nfevals,...
                        npds, ndecomps, nsolves],...
                        nout, tout, yout,...
                        haveEventFcn, teout, yeout, ieout,...
                        {t2data,y2data,idxNonNegative},...
                        t);
                    if nargout > 0
                        varargout = solver_output;
                    end
                    return;
                end

                nofailed = false;
                abshlast = absh;
                if NNrejectStep
                    absh = max(hmin, 0.5*absh);
                else
                    absh = max(abshlast * max(0.1, 0.7*(rtol/err)^pow), hmin);
                end
                h = tdir * absh;
                z = (absh / abshlast) * z;
                needNewLU = true;
                done = false;
            else                              % Successful step
                break;

            end
        end
    end % while true
    nsteps = nsteps + 1;

    NNreset_znew = false;
    if nonNegative && any(ynew(idxNonNegative) < 0)
        NNidx = idxNonNegative(ynew(idxNonNegative) < 0); % logical indexing
        ynew(NNidx) = 0;
        if normcontrol
            normynew = norm(ynew);
        end
        NNreset_znew = true;
    end

    if haveEventFcn
        [te,ye,ie,valt,stop] = odezero(@ntrp23tb,eventFcn,eventArgs,valt,...
            t,y,tnew,ynew,t0,t2,y2,idxNonNegative);

        if ~isempty(te)
            if output_sol || (nargout > 2)
                teout = [teout, te]; %#ok<AGROW>
                yeout = [yeout, ye]; %#ok<AGROW>
                ieout = [ieout, ie]; %#ok<AGROW>
            end
            if stop               % Stop on a terminal event.
                % Adjust the interpolation data to [t te(end)].
                t2_zc = t + alpha*(te(end)-t);
                y2_zc = ntrp23tb(t2_zc,t,y,tnew,ynew,t2,y2,idxNonNegative);
                t2 = t2_zc;
                y2 = y2_zc;
                tnew = te(end);
                ynew = ye(:,end);
                done = true;
            end
        end
    end

    if output_sol
        nout = nout + 1;
        if nout > length(tout)
            tout = [tout, zeros(1,chunk)]; %#ok<AGROW> requires chunk >= refine
            yout = [yout, zeros(neq,chunk)]; %#ok<AGROW>
            t2data = [t2data, zeros(1,chunk)]; %#ok<AGROW>
            y2data = [y2data, zeros(neq,chunk)]; %#ok<AGROW>
        end
        tout(nout) = tnew; %#ok<AGROW>
        yout(:,nout) = ynew; %#ok<AGROW>
        t2data(nout) = t2; %#ok<AGROW>
        y2data(:,nout) = y2; %#ok<AGROW>
    end

    if output_ty || haveOutputFcn
        switch outputAt
            case 'SolverSteps'        % computed points, no refinement
                nout_new = 1;
                tout_new = tnew;
                yout_new = ynew;
            case 'RefinedSteps'       % computed points, with refinement
                tref = t + (tnew-t)*S;
                nout_new = refine;
                tout_new = [tref, tnew];
                yout_new = [ntrp23tb(tref,t,y,tnew,ynew,t2,y2,idxNonNegative), ynew];
            case 'RequestedPoints'    % output only at tspan points
                nout_new =  0;
                tout_new = [];
                yout_new = [];
                while next <= ntspan
                    if tdir * (tnew - tspan(next)) < 0
                        if haveEventFcn && stop     % output tstop,ystop
                            nout_new = nout_new + 1;
                            tout_new = [tout_new, tnew]; %#ok<AGROW>
                            yout_new = [yout_new, ynew]; %#ok<AGROW>
                        end
                        break;
                    end
                    nout_new = nout_new + 1;
                    tout_new = [tout_new, tspan(next)]; %#ok<AGROW>
                    if tspan(next) == tnew
                        yout_new = [yout_new, ynew]; %#ok<AGROW>
                    else
                        yout_new = [yout_new, ntrp23tb(tspan(next),t,y,tnew,ynew,t2,y2,idxNonNegative)]; %#ok<AGROW>
                    end
                    next = next + 1;
                end
        end

        if nout_new > 0
            if output_ty
                oldnout = nout;
                nout = nout + nout_new;
                if nout > length(tout)
                    tout = [tout, zeros(1,chunk)]; %#ok<AGROW> requires chunk >= refine
                    yout = [yout, zeros(neq,chunk)]; %#ok<AGROW>
                end
                idx = oldnout+1:nout;
                tout(idx) = tout_new; %#ok<AGROW>
                yout(:,idx) = yout_new; %#ok<AGROW>
            end
            if haveOutputFcn
                stop = feval(outputFcn,tout_new,yout_new(outputs,:),'',outputArgs{:});
                if stop
                    done = true;
                end
            end
        end
    end

    if done
        break
    end

    % Advance the integration one step.
    if NNreset_znew
        % Used znew for unperturbed solution to interpolate.  In perturbing ynew,
        % defined NNidx.  Use now to reset znew to move along constraint.
        znew(NNidx) = 0;
    end
    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,4} % state dependent
            % 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 % while ~done

solver_output = odefinalize(solver_name, sol,...
    outputFcn, outputArgs,...
    printstats, [nsteps, nfailed, nfevals,...
    npds, ndecomps, nsolves],...
    nout, tout, yout,...
    haveEventFcn, teout, yeout, ieout,...
    {t2data,y2data,idxNonNegative});
if nargout > 0
    varargout = solver_output;
end

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

function [y,z,iter,itfail,rate] = ...
    itsolve(M,t,y,z,d,h,Factors,piv,ode,rtol,wt,rate,Mtype,Mfun,Margs,odeArgs)
% Solve the nonlinear equation M*z = h*f(t,v+d*z) and y = v+d*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 given by ode(t,y,odeArgs{:}).
% Similarly, if M is state-dependent, it is given by Mfun(t,y,Margs{:}).
% Factors,piv is the lu decomposition of the matrix M-d*h*dfdy, where dfdy is an
% approximate Jacobian of f. For full matrices, P is the permutation vector
% and Q and R are left empty. 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  % state dependent
        M = Mfun(t,y,Margs{:});
    end
    rhs = h * ode(t,y,odeArgs{:}) - M * z;

    del = odesolve(Factors, piv, rhs);

    z = z + del;
    y = y + d*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

