function [ tvals, yvals ] = hamming ( f, tspan, startval, nstep )
%
%  function [ tvals, yvals ] = hamming ( f, tspan, startval, nstep )
%
%  HAMMING uses Hamming's method to solve an ODE.
%
%  F is the name of the user-supplied right hand side function.
%  TSPAN(1) is the starting time, TSPAN(2) the stopping time.
%  STARTVAL is the column vector of initial values for the solution.
%  NSTEP is the number of steps to take from start to finish.
%
%  Reference:
%    George Lindfield and John Penny,
%    Numerical Methods using MATLAB, second edition,
%    Prentice Hall, 2000, pages 219-220.
%
b = [];
c = [];
d = [];
order = 4;

b = [ 1/6, 1/3, 1/3, 1/6];
c = [ 0.0, 0.0, 0.0, 0.0;
      0.5, 0.0, 0.0, 0.0;
      0.0, 0.5, 0.0, 0.0;
      0.0, 0.0, 1.0, 0.0 ];
d = [0.0, 0.5, 0.5, 1.0];

h = ( tspan(2) - tspan(1) ) / nstep;
y = startval;
t = tspan(1);
ys(1) = startval;

yvals = startval;
tvals = tspan(1);
%
%  Use a Runge Kutta method to get the early values.
%
for j = 2 : 4
  k(1) = h * feval ( f, t, y );
  for i = 2 : order
    k(i) = h * feval ( f, t + h * d(i), y + c(i,1:i-1) * k(1:i-1)' );
  end;
  y1 = y + b * k';
  ys(j) = y1;
  t1 = t + h;
  fval(j) = feval ( f, t1, y1 );
  tvals = [ tvals, t1 ];
  yvals = [ yvals, y1 ];
  t = t1;
  y = y1;
end;
%
%  Once you have values for at least 4 times steps, switch to 
%  Hamming's method.
%
%  PLEASE FIX THIS.  C IS USED TWO PLACES WITH TWO MEANINGS.
%
for i = 5 : nstep+1

  y1 = ys(1) + 4.0 * h * ( 2.0 * fval(4) - fval(3) + 2.0 * fval(2) ) / 3.0;
  t1 = t + h;

  if ( i = 5 ) 
    y1m = y1;
  else
    y1m = y1 + 112.0 * ( c - p ) / 121.0;
  end;

  fval(5) = feval ( f, t1, y1m );
  yc = ( 9.0 * ys(4) - ys(2) + 3.0 * h * ( 2.0 * fval(4) + fval(5) - fval(3) ) ) / 8.0;
  ycm = yc + 9.0 * ( y1 - yc ) / 121.0;
  p = y1;
  c = yc;
%
%  Isn't this the damndest goofy update list you ever saw?
%
  fval(5) = feval ( f, t1, ycm );
  fval(2) = fval(3);
  fval(3) = fval(4);
  fval(4) = fval(5);

  ys(1) = ys(2);
  ys(2) = ys(3);
  ys(3) = ys(4);
  ys(4) = ycm;
%
%  Add T1, YCM to the solution vector.
%
  tvals = [ tvals, t1 ];
  yvals = [ yvals, ycm ];
  t = t1;
end;
