function [tout,uout] = ode23smp(f,tspan,u0,tol)
% Simplified version of Matlab's ode23
% Integrates the ODE system du/dt = f(t,u) from t=t0 to t=T
%   with initial condition u(t0)=u0 
%   using a 3-stage Runge-Kutta method (due to Bogacki and Shampine)
%   with adaptive step size h
% If the 4th input argument tol is a number, this is taken to be
%   the relative error tolerance (in place of the default value 1e-3);
% if the 4th input argument is of the form 
%   tol=odeset('RelTol',RTOL,'AbsTol',ATOL) where RTOL,ATOL are numbers,
%   then the program uses RTOL as the relative error tolerance
%   and ATOL as the absolute error tolerance (instead of the default 1e-6)

% default tolerances
rtol = 1.e-3;
atol = 1.e-6;
if nargin == 4 && isnumeric(tol)
   rtol = tol;
elseif nargin == 4 && isstruct(tol)
   if ~isempty(tol.RelTol), rtol = tol.RelTol; end
   if ~isempty(tol.AbsTol), atol = tol.AbsTol; end
end
threshold = atol/rtol;

% initialize
t0 = tspan(1);
T = tspan(2);
hmax = 0.1*(T-t0);
t = t0;
u = u0(:);
tout(1)=t;
uout(1,:)=u.';

% choose an initial step size based on the overall scale of the problem
s1 = f(t,u); % calculate the first slope
r = norm(s1./max(abs(u),threshold), inf) + realmin;
h = 0.8*rtol^(1/3)/r;

% the time-stepping loop
while t < T  
   hmin = 16*eps*abs(t);
   if h > hmax, h = hmax; end
   if h < hmin, h = hmin; end   
   if 1.1*h >= T - t % stretch the step if t is close to T
      h = T - t;
   end
   % take a step (s1 was already calculated above)
   s2 = f(t+h/2, u+h/2*s1);
   s3 = f(t+3*h/4, u+3*h/4*s2);
   tnew = t + h;
   unew = u + h*(2*s1 + 3*s2 + 4*s3)/9;
   % estimate the error
   s4 = f(tnew, unew);
   e = h*(-5*s1 + 6*s2 + 8*s3 - 9*s4)/72;
   err_est = norm(e./max( max(abs(u),abs(unew)), threshold ), inf) + realmin;
   % accept the solution if the estimated error is less than the tolerance
   if err_est <= rtol
      t = tnew;
      u = unew; 
      tout(end+1,1) = t;
      uout(end+1,:) = u.';
      s1 = s4; 
   end
   % find a new step size
   h = h*min(5, 0.8*(rtol/err_est)^(1/3)); 
   % exit early if step size is too small
   if h <= hmin
      warning('Step size %e too small at t = %e.\n',h,t);
      break
   end
end

end
