function [tout,yout,stats] = ode45(ydot,tspan,y0,options)
%ODE45	Solve non-stiff differential equations, medium order method.
%	[T,Y] = ODE45('ydot',TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates
%	the system of first order differential equations y' = ydot(t,y) from
%	time T0 to TFINAL with initial conditions Y0.  Function ydot(t,y)
%	must return a column vector.  Each row in solution matrix Y
%	corresponds to a time returned in column vector T.  To obtain
%	solutions at the specific times T0, T1, ..., TFINAL (all increasing
%	or all decreasing), use TSPAN = [T0 T1 ... TFINAL].
%	
%	[T,Y] = ODE45('ydot',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 'rtol' (1e-3 by
%	default) and vector of absolute error tolerances 'atol' (all
%	components 1e-6 by default).
%	
%	It is possible to specify tspan, y0 and options in ydot.  If TSPAN
%	or Y0 is empty, or if ODE45 is invoked as ODE45('ydot'), ODE45 calls
%	[tspan,y0,options] = ydot([],[]) to obtain any values not supplied
%	at the command line.  TYPE R3BODY to see how this is coded.
%	
%	As an example, the commands
%	
%	    options = odeset('rtol',1e-4,'atol',[1e-4 1e-4 1e-5]);
%	    ode45('rigid',[0 12],[0 1 1],options);
%	
%	solve the system y' = rigid(t,y) with relative error tolerance 1e-4
%	and absolute tolerances of 1e-4 for the first two components and
%	1e-5 for the third.  When called with no output arguments, as in
%	this example, ODE45 calls the default output function ODEPLOT to
%	plot the solution as it is computed.
%	
%	See also
%	    other ODE solvers:   ODE23, ODE113, ODE15S, ODE23S
%	    options handling:    ODESET, ODEGET
%	    output functions:    ODEPLOT, ODEPHAS2, ODEPHAS3
%	    ydot examples:       R3BODY, RIGID, TWOBODY, VDPNS

%	ODE45 is an implementation of the explicit Runge-Kutta (4,5) pair of
%	Dormand and Prince called variously RK5(4)7FM, DOPRI5, DP(4,5) and
%	DP54.  It uses a "free" interpolant of order 4 communicated
%	privately by Dormand and Prince.  Local extrapolation is done.
%	Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
%	M. W. Reichelt, Rept. 94-6, Math. Dept., SMU, Dallas, TX, 1994.

%	Mark W. Reichelt and Lawrence F. Shampine, 6-14-94
%	Copyright (c) 1984-95 by The MathWorks, Inc.

false = 0;
true = ~false;

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

if nargin == 1
  tspan = []; y0 = []; options = [];
elseif nargin == 2
  y0 = []; options = [];
elseif nargin == 3
  options = [];
end

% Get default tspan and y0 from ydot if none are specified.
if (length(tspan) == 0) | (length(y0) == 0)
  [ydot_tspan,ydot_y0,ydot_options] = feval(ydot,[],[]);
  if length(tspan) == 0
    tspan = ydot_tspan;
  end
  if length(y0) == 0
    y0 = ydot_y0;
  end
  if length(options) == 0
    options = ydot_options;
  else
    options = odeset(ydot_options,options);
  end
end

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

y = y0(:);
neq = length(y);

% Get options, and set defaults.
rtol = odeget(options,'rtol',1e-3);
atol = odeget(options,'atol',1e-6);
atol = atol(:);
if (rtol <= 0) | any(atol <= 0)  
  error('Tolerance rtol and all entries of atol must be positive.');
end
if rtol < 100 * eps 
  rtol = 100 * eps;
  fprintf('Warning: rtol has been increased to %e\n',rtol);
end
if length(atol) == 1
  atol = atol + zeros(neq,1);
elseif length(atol) ~= neq
  error(['Vector atol must be same length as y0 vector (' num2str(neq) ').']);
end
threshold = atol ./ rtol;

userhmax = abs(odeget(options,'hmax',0.1*abs(tfinal-t)));
hmax = min(userhmax, abs(tfinal-t));

hmin = 40 * eps * abs(t);
userhmin = abs(odeget(options,'hmin',hmin));
hmin = min(max(hmin, userhmin), hmax);

if nargout ~= 0
  outfun = odeget(options,'outfun');
else
  outfun = odeget(options,'outfun','odeplot');
end

refine = odeget(options,'refine',4);    % Necessary for smooth plots.

printstats = odeget(options,'printstats',false);

if length(odeget(options,'mass')) ~= 0
  error('Solver does not handle mass matrices, M*y''.  See ODE15S or ODE23S.');
end

% Initialize the output function.
if length(outfun) == 0
  haveoutfun = false;
else
  haveoutfun = true;
  feval(outfun,[t tfinal],y,'init');
end

% Allocate memory if we're generating output.
if nargout ~= 0
  if (ntspan == 1) & (refine == 0) 	% only 1 output at tfinal
    nout = 0;
  else
    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
end

% Initialize method parameters.
pow = 1/5;
A = [1/5; 3/10; 4/5; 8/9; 1; 1];
B = [
    1/5		3/40	44/45	19372/6561	9017/3168	35/384
    0		9/40	-56/15	-25360/2187	-355/33		0
    0		0	32/9	64448/6561	46732/5247	500/1113
    0		0	0	-212/729	49/176		125/192
    0		0	0	0		-5103/18656	-2187/6784
    0		0	0	0		0		11/84
    0		0	0	0		0		0
    ];
E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];
BI = [
    1		-183/64		37/12		-145/128
    0		0		0		0
    0		1500/371	-1000/159	1000/371
    0		-125/32		125/12		-375/64	
    0		9477/3392	-729/106	25515/6784
    0		-11/7		11/3		-55/28
    0		3/2		-4		2.5
    ];
f = zeros(neq,7);

% Set the output flag.
if (ntspan > 2) | (refine == 0)
  outflag = 1; 				% output only at tspan points
  onesS = ones(4,1);
elseif refine == 1
  outflag = 2; 				% computed points, no refinement
else
  outflag = 3; 				% computed points, with refinement
  ones1r = ones(1,refine-1);
  S = (1:refine-1) / refine;
  S = cumprod(S(ones(4,1),:));
end

% Compute an initial step size h using y'(t).
f0 = feval(ydot,t,y);
[m,n] = size(f0);
if n > 1
  error(['Function ' ydot '(t,y) must return a column vector.'])
elseif m ~= neq
  error(['Vector ' ydot '(t,y) must be same length as initial conditions.']);
end
nfevals = nfevals + 1; 			% stats
absh = min(hmax, abs(tspan(next) - t));
rh = norm(f0 ./ (abs(y) + threshold),inf) / (0.8 * rtol^pow);
if absh * rh > 1
  absh = 1 / rh;
end
absh = max(absh, hmin);
h = tdir * absh;
f(:,1) = f0;

% THE MAIN LOOP

done = false;
while ~done
  
  hmin = max(userhmin, 40 * eps * abs(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
    
    hA = h * A;
    hB = h * B;
    f(:,2) = feval(ydot, t + hA(1), y + f*hB(:,1));
    f(:,3) = feval(ydot, t + hA(2), y + f*hB(:,2));
    f(:,4) = feval(ydot, t + hA(3), y + f*hB(:,3));
    f(:,5) = feval(ydot, t + hA(4), y + f*hB(:,4));
    f(:,6) = feval(ydot, t + hA(5), y + f*hB(:,5));
    tnew = t + hA(6);
    ynew = y + f*hB(:,6);
    f(:,7) = feval(ydot, tnew, ynew);
    nfevals = nfevals + 6; 		% stats
    
    % Estimate the error.
    err = absh * norm((f * E) ./ (max(abs(y),abs(ynew)) + threshold),inf);
    
    % 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
      nfailed = nfailed + 1; 		% stats
      if absh <= hmin
	fprintf('Step failure at %e with a minimum step size of %e\n', t, h);
	if nargout ~= 0
	  tout = tout(1:nout);
	  yout = yout(1:nout,:);
	  stats = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
	end
	return;
      end
      
      if nofailed
	nofailed = false;
	absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
      else
        absh = max(hmin, 0.5 * absh);
      end
      h = tdir * absh;
      done = false;
      
    else 				% Successful step
      break;
      
    end
  end
  
  if nargout ~= 0
    oldnout = nout;
    if 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
      hS = h * S;
      i = oldnout+1:nout-1;
      tout(i) = t + hS(1,:)';
      yout(i,:) = (y(:,ones1r) + f*BI*hS).';
      tout(nout) = tnew;
      yout(nout,:) = ynew.';
    elseif 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 == 1 		% output only at tspan points
      c = f * (h * BI);
      while true
	if next > ntspan 
	  break;
	elseif tdir * (tnew - tspan(next)) < 0
	  break;
	end
	nout = nout + 1; 		% tout and yout are already allocated
	tout(nout) = tspan(next);
	s = (tspan(next) - t) / h;
	yout(nout,:) = (y + c*cumprod(s(onesS))).';
	next = next + 1;
      end
    end
    
    if haveoutfun
      for i = oldnout+1:nout
	if feval(outfun,tout(i),yout(i,:).') == 1
	  tout = tout(1:nout);
	  yout = yout(1:nout,:);
	  stats = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
	  return;
	end
      end
    end
    
  elseif haveoutfun
    if outflag == 3 			% computed points, with refinement
      hS = h * S;
      for i = 1:refine-1
	yinterp = y + f*BI*hS(:,i);
	if feval(outfun,t + hS(1,i),yinterp) == 1
	  return;
	end
      end
      if feval(outfun,tnew,ynew) == 1
	return;
      end
    elseif outflag == 2
      if feval(outfun,tnew,ynew) == 1
	return;
      end
    elseif outflag == 1 		% output only at tspan points
      c = f * (h * BI);
      while true
	if next > ntspan 
	  break;
	elseif tdir * (tnew - tspan(next)) < 0
	  break;
	end
	s = (tspan(next) - t) / h;
	yinterp = y + c*cumprod(s(onesS));
	if feval(outfun,tspan(next),yinterp) == 1
	  return;
	end
	next = next + 1;
      end
    end
  end
  
  % If there were no failures compute a new h.
  if nofailed 				% note h may shrink by 0.8
    absh = min(hmax, absh / max(0.2, 1.25*(err/rtol)^pow)); % err may be 0
    h = tdir * absh;
  end
  
  % Advance the integration one step.
  t = tnew;
  y = ynew;
  nsteps = nsteps + 1; 			% stats
  f(:,1) = f(:,7); 			% Already evaluated ydot at (tnew,ynew)
  
end

if printstats 				% print cost statistics
  fprintf('%g successful steps\n', nsteps);
  fprintf('%g failed attempts\n', nfailed);
  fprintf('%g calls to ydot\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,:);
  stats = [nsteps; nfailed; nfevals; npds; ndecomps; nsolves];
end
