function [t,u,e] = ssode23(neqn,t0,u0,h)
%
% Function ssode23 computes an ODE solution by the embedded
% RK (2,3) pair discussed by Shampine, L. F., I. Gladwell
% and S. Thompson, Solving ODEs with Matlab, Cambridge
% University Press, Cambridge, New York, 2003, p 54, Table
% 2.3.  It also estimates the truncation error of the
% solution, and applies this estimate as a correction to
% the solution vector.
%
% Argument list
%
%   neqn     number of first order ODEs
%
%   t0       initial value of independent variable
%
%   u0       initial condition vector of length neqn
%
%   h        integration step
%
%   t        independent variable
%
%   u        ODE solution vector of length neqn after
%            one ode23 step
%   
%   e        estimate of truncation error of the solu-
%            tion vector
%
% The following information is taken from the source code
% of the Matlab integrator ode23:
%
% Initialize method parameters.
% 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];
%
% This information can be compared with the coding below
% of the RK (2,3) pair.
%
% Derivative vector at initial (base) point
  [ut0]=derv(neqn,t0,u0);
%
% k1, advance of dependent variable vector and
% independent variable for calculation of k2
  k1=h*ut0;
  u=u0+0.5*k1;
  t=t0+0.5*h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k2, advance of dependent variable vector and
% independent variable for calculation of k3
  k2=h*ut;
  u=u0+0.75*k2;
  t=t0+0.75*h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k3, advance of dependent variable vector and
% independent variable for calculation of k4
  k3=h*ut;
  u=u0+2.0/9.0*k1+1.0/3.0*k2+4.0/9.0*k3;
  t=t0+h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k4
  k4=h*ut;
%
% Second order step
  sum2=u0+7.0/24.0*k1+1.0/4.0*k2+1.0/3.0*k3+1.0/8.0*k4;
%
% Third order step
  sum3=u;
  t=t0+h;
%
% Truncation error estimate
% e = sum3 - sum2
%
%   = u0+2.0/9.0*k1+1.0/3.0*k2+4.0/9.0*k3
%
%   - u0+7.0/24.0*k1+1.0/4.0*k2+1.0/3.0*k3+1.0/8.0*k4
%
%   = -5.0/72.0*k1+1.0/12.0*k2+1.0/9.0*k3-1.0/8.0*k4
%
% which reflects the weighting coefficients in the vector
% E indicated above.
  e=sum3-sum2;
%
% End of ssode23
