function [t,u,e] = ssode45(neqn,t0,u0,h)
%
% Function ssode45 computes an ODE solution by the embedded
% RK (4,5) pair discussed by Shampine, L. F., I. Gladwell
% and S. Thompson, Solving ODEs with Matlab, Cambridge
% University Press, Cambridge, New York, 2003, p 53.  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 ode45 step
%   
%   e        estimate of truncation error of the solu-
%            tion vector
%
% The following information is taken from the source code
% of the Matlab integrator ode45:
%
% 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];
%
% This information can be compared with the coding below
% of the RK (4,5) 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+(1.0/5.0)*k1;
  t=t0+(1.0/5.0)*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+(3.0/40.0)*k1...
      +(9.0/40.0)*k2;
  t=t0+(3.0/10.0)*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+(44.0/45.0)*k1...
      -(56.0/15.0)*k2...
      + (32.0/9.0)*k3;
  t=t0+(4.0/5.0)*h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k4, advance of dependent variable vector and
% independent variable for calculation of k5
  k4=h*ut;
  u=u0+(19372.0/6561.0)*k1...
      -(25360.0/2187.0)*k2...
      +(64448.0/6561.0)*k3...
      -   (212.0/729.0)*k4;
  t=t0+(8.0/9.0)*h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k5, advance of dependent variable vector and
% independent variable for calculation of k6
  k5=h*ut;
  u=u0+ (9017.0/3168.0)*k1...
      -    (355.0/33.0)*k2...
      +(46732.0/5247.0)*k3...
      +    (49.0/176.0)*k4...
      -(5103.0/18656.0)*k5;
  t=t0+h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k6, advance of dependent variable vector and
% independent variable for calculation of k7
% (u is the fifth-order result)
  k6=h*ut;
  u=u0+   (35.0/384.0)*k1...
      + (500.0/1113.0)*k3...
      +  (125.0/192.0)*k4...
      -(2187.0/6784.0)*k5...
      +    (11.0/84.0)*k6;
  t=t0+h;
%
% Derivative vector at new u, t
  [ut]=derv(neqn,t,u);
%
% k7
  k7=h*ut;
%
% Truncation error estimate for u
  e=     (71.0/57600.0)*k1...
    -    (71.0/16695.0)*k3...
    +     (71.0/1920.0)*k4...
    -(17253.0/339200.0)*k5...
    +      (22.0/525.0)*k6...
    -        (1.0/40.0)*k7;
%
% End of ssode45
