function [tout,yout,varargout] = ode23(ode,tspan,y0,options,varargin)
%ODE23  Solve non-stiff differential equations, low order method.
%   [T,Y] = ODE23(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. Function ODEFUN(T,Y) must return a column vector
%   corresponding to f(t,y). Each row in the solution array Y corresponds to
%   a time returned in the column vector T. To obtain solutions at specific
%   times T0,T1,...,TFINAL (all increasing or all decreasing), use 
%   TSPAN = [T0 T1 ... TFINAL].   
%   
%   [T,Y] = ODE23(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
%   integration properties 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] = ODE23(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2,...) passes the additional
%   parameters P1,P2,... to the ODE function as ODEFUN(T,Y,P1,P2...), and to
%   all functions specified in OPTIONS. Use OPTIONS = [] as a place holder if
%   no options are set.   
%   
%   ODE23 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
%   nonsingular. Use ODESET to set the 'Mass' property to a function 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' property. 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'. ODE15S and ODE23T can solve problems with singular mass matrices.  
%
%   [T,Y,TE,YE,IE] = ODE23(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
%   property in OPTIONS set to a function 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 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.   
%
%   Example    
%         [t,y]=ode23(@vdp1,[0 20],[2 0]);   
%         plot(t,y(:,1));
%     solves the system y' = vdp1(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 
%       other ODE solvers:  ODE45, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB
%       options handling:   ODESET, ODEGET
%       output functions:   ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT
%       ODE examples:       RIGIDODE, BALLODE, ORBITODE
%
%   NOTE: 
%     The interpretation of the first input argument of the ODE solvers and
%     some properties available through ODESET have changed in this version
%     of MATLAB. Although we still support the v5 syntax, any new
%     functionality is available only with the new syntax. To see the v5
%     help, type in the command line 
%         more on, type ode23, more off

%   NOTE:
%     This portion describes the v5 syntax of ODE45.
%
%   [T,Y] = ODE23('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] = ODE23('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] = ODE23('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 ODE23 calls the ODE file
%   [TSPAN,Y0,OPTIONS] = F([],[],'init') to obtain any values not supplied
%   in the ODE23 argument list.  Empty arguments at the end of the call list
%   may be omitted, e.g. ODE23('F').
%   
%   ODE23 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.
%   ODE15S and ODE23T can solve problems with singular mass matrices.
%   
%   [T,Y,TE,YE,IE] = ODE23('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 

%   ODE23 is an implementation of the explicit Runge-Kutta (2,3) pair of
%   Bogacki and Shampine called BS23.  It uses a "free" interpolant of order
%   3.  Local extrapolation is done.

%   Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
%   M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997.

%   Mark W. Reichelt and Lawrence F. Shampine, 6-14-94
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 5.68 $  $Date: 2000/08/01 20:57:29 $

true = 1;
false = ~true;

stats = struct('nsteps',0,'nfailed',0,'nfevals',0,...
               'npds',0,'ndecomps',0,'nsolves',0); 

if nargin < 4
  options = [];
  if nargin < 3
    y0 = [];
    if nargin < 2
      tspan = [];
      if nargin < 1
        error('Not enough input arguments.  See ODE23.');
      end  
    end
  end
end
  
FcnHandlesUsed = isa(ode,'function_handle');

% Handle solver arguments
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, ...
          options, atol, rtol, threshold, normcontrol, normy, hmax, htry]  ...
         = odearguments(FcnHandlesUsed,'ode23', ode, tspan, y0, options, ...
                        varargin);
stats.nfevals = stats.nfevals + 1;

% Handle the output
if nargout > 0
  outfun = odeget(options,'OutputFcn',[],'fast');
else
  outfun = odeget(options,'OutputFcn',@odeplot,'fast');
end
if isempty(outfun)
  haveoutfun = false;
else
  haveoutfun = true;
  outputs = odeget(options,'OutputSel',1:neq,'fast');
  if isa(outfun,'function_handle')
    outputArgs = [{''},varargin];
    outputArgs1 = varargin;
  else       % With v5 syntax do not pass additional input arguments to outfun
    outputArgs = {};      
    outputArgs1 = {};
  end  
end
refine = odeget(options,'Refine',1,'fast');
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
printstats = strcmp(odeget(options,'Stats','off','fast'),'on');

% Handle the event function 
[haveeventfun,eventfun,eventargs,valt,teout,yeout,ieout] = ...
    odeevents(FcnHandlesUsed,ode,t0,y0,options,varargin);

% Handle the mass matrix 
[Mtype, Mfun, Margs, M] =  odemass(FcnHandlesUsed,ode,t0,y0,options,varargin);
if Mtype>0  % non-trivial mass matrix
  Msingular = odeget(options,'MassSingular','no','fast');
  if strcmp(Msingular,'maybe')
    warning(['This solver assumes MassSingular is ''no''.  See ODE15S or ' ...
             'ODE23T.']);
  elseif strcmp(Msingular,'yes')
     error(['MassSingular cannot be ''yes'' for this solver.  See ODE15S '...
            ' or ODE23T.']);
  end
  if Mtype == 1
    [L,U] = lu(M);
  else
    L = [];
    U = [];
  end  
  [ode,args] = odemassexplicit(FcnHandlesUsed,Mtype,ode,args,Mfun,Margs,L,U);

  % Recompute f0, this time fun takes care about the mass matrix
  f0 = feval(ode,t0,y0,args{:});
  stats.nfevals = stats.nfevals + 1;
end

t = t0;
y = y0;

% 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 = min(max(100,50*refine),floor((2^13)/neq));
    tout = zeros(chunk,1);
    yout = zeros(chunk,neq);
  end
  nout = 1;
  tout(nout) = t;
  yout(nout,:) = y.';
end

% Initialize method parameters.
pow = 1/3;
A = [1/2; 3/4; 1];
B = [
    1/2         0               2/9
    0           3/4             1/3
    0           0               4/9
    0           0               0
    ]; 
E = [-5/72; 1/12; 1/9; -1/8];
f = zeros(neq,4);
hmin = 16*eps*abs(t);
if isempty(htry)
  % Compute an initial step size h using y'(t).
  absh = min(hmax, abs(tspan(next) - t));
  if normcontrol
    rh = (norm(f0) / max(normy,threshold)) / (0.8 * rtol^pow);
  else
    rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
  end
  if absh * rh > 1
    absh = 1 / rh;
  end
  absh = max(absh, hmin);
else
  absh = min(hmax, max(hmin, htry));
end
f(:,1) = f0;

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

% THE MAIN LOOP

done = false;
while ~done
  
  % By default, hmin is a small number such that t+hmin is only slightly
  % different than t.  It might be 0 if t is 0.
  hmin = 16*eps*abs(t);
  absh = min(hmax, max(hmin, absh));    % couldn't limit absh until new hmin
  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
  
  % LOOP FOR ADVANCING ONE STEP.
  nofailed = true;                      % no failed attempts
  while true
    hB = h * B;
    f(:,2) = feval(ode,t+h*A(1),y+f*hB(:,1),args{:});
    f(:,3) = feval(ode,t+h*A(2),y+f*hB(:,2),args{:});
    tnew = t + h*A(3);
    if done
      tnew = tfinal;   % Hit end point exactly.
    end
    ynew = y + f*hB(:,3);
    f(:,4) = feval(ode,tnew,ynew,args{:});
    stats.nfevals = stats.nfevals + 3;              
    
    % Estimate the error.
    if normcontrol
      normynew = norm(ynew);
      err = absh * (norm(f * E) / max(max(normy,normynew),threshold));
    else
      err = absh * norm((f * E) ./ max(max(abs(y),abs(ynew)),threshold),inf);
    end
    
    % Accept the solution only if the weighted error is no more than the
    % tolerance rtol.  Estimate an h that will yield an error of rtol on
    % the next step or the next try at taking this step, as the case may be,
    % and use 0.8 of this value to avoid failures.
    if err > rtol                       % Failed step
      stats.nfailed = stats.nfailed + 1;            
      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',outputArgs1{:});
        end
        if printstats                   % print cost statistics
          fprintf('%g successful steps\n', stats.nsteps);
          fprintf('%g failed attempts\n', stats.nfailed);
          fprintf('%g function evaluations\n', stats.nfevals);
          fprintf('%g partial derivatives\n', stats.npds);
          fprintf('%g LU decompositions\n', stats.ndecomps);
          fprintf('%g solutions of linear systems\n', stats.nsolves);
        end
        if nargout > 0
          tout = tout(1:nout);
          yout = yout(1:nout,:);
          if haveeventfun
            varargout{1} = teout;
            varargout{2} = yeout;
            varargout{3} = ieout;
            if ~FcnHandlesUsed
              varargout{4} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                              stats.npds; stats.ndecomps; stats.nsolves];
            end
          else
            if ~FcnHandlesUsed
              varargout{1} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                              stats.npds; stats.ndecomps; stats.nsolves];
            end          
          end
        end
        return;
      end
      
      if nofailed
        nofailed = false;
        absh = max(hmin, absh * max(0.5, 0.8*(rtol/err)^pow));
      else
        absh = max(hmin, 0.5 * absh);
      end
      h = tdir * absh;
      done = false;
      
    else                                % Successful step
      break;
      
    end
  end
  stats.nsteps = stats.nsteps + 1;                  
  
  if haveeventfun
    [te,ye,ie,valt,stop] = ...
        odezero(@ntrp23,eventfun,eventargs,valt,t,y,tnew,ynew,t0,h,f);
    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)];  % requires chunk >= refine
        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,:) = ntrp23(tout(i),t,y,[],[],h,f).';
      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,:) = ntrp23(tspan(next),t,y,[],[],h,f).';
        next = next + 1;
      end
    end
    
    if haveoutfun
      i = oldnout+1:nout;
      if ~isempty(i) & (feval(outfun,tout(i),yout(i,outputs).',outputArgs{:}) == 1)
        tout = tout(1:nout);
        yout = yout(1:nout,:);
        if haveeventfun
          varargout{1} = teout;
          varargout{2} = yeout;
          varargout{3} = ieout;
          if ~FcnHandlesUsed
            varargout{4} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                            stats.npds; stats.ndecomps; stats.nsolves];
          end
        else
          if ~FcnHandlesUsed
            varargout{1} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                            stats.npds; stats.ndecomps; stats.nsolves];
          end       
        end
        return;
      end
    end
    
  elseif haveoutfun
    if outflag == 2
      if feval(outfun,tnew,ynew(outputs),outputArgs{:}) == 1
        return;
      end
    elseif outflag == 3                 % computed points, with refinement
      tinterp = t + (tnew-t)*S;
      yinterp = ntrp23(tinterp,t,y,[],[],h,f);
      if feval(outfun,[tinterp; tnew],[yinterp(outputs,:), ynew(outputs)],outputArgs{:}) == 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) = ntrp23(tspan(next),t,y,[],[],h,f);
        next = next + 1;
      end
      if ninterp > 0
        if feval(outfun,tinterp(1:ninterp),yinterp(outputs,1:ninterp),outputArgs{:}) == 1
          return;
        end
      end
    end
  end
  
  % If there were no failures compute a new h.
  if nofailed
    % Note that absh may shrink by 0.8, and that err may be 0.
    temp = 1.25*(err/rtol)^pow;
    if temp > 0.2
      absh = absh / temp;
    else
      absh = 5.0*absh;
    end
  end
  
  % Advance the integration one step.
  t = tnew;
  y = ynew;
  if normcontrol
    normy = normynew;
  end
  f(:,1) = f(:,4);                      % Already evaluated ode(tnew,ynew)
  
end

if haveoutfun
  feval(outfun,[],[],'done',outputArgs1{:});
end

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

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

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

function yinterp = ntrp23(tinterp,t,y,tnew,ynew,h,f)
%NTRP23 Interpolation helper function for ODE23.
%   YINTERP = NTRP23(TINTERP,T,Y,TNEW,YNEW,H,F) uses data computed in ODE23
%   to approximate the solution at time TINTERP.

BI = [
    1       -4/3        5/9
    0       1       -2/3
    0       4/3     -8/9
    0       -1      1
    ];

s = ((tinterp - t) / h)';       % may be a row vector

yinterp = y(:,ones(length(tinterp),1)) + f*(h*BI)*cumprod(s(ones(3,1),:));
