function [ x, y ] = ab2 ( f, x_range, y_initial, nstep )
%
%  function [ x, y ] = ab2 ( f, x_range, y_initial, nstep )
%
%  AB2 uses NSTEP steps of the Adams-Bashforth method of order 2 to 
%  estimate Y, the solution of an ODE, at the equally spaced points X in 
%  the range X_RANGE(1) to X_RANGE(2).  The name of the derivative function 
%  is F. 
%
dx = ( x_range(2) - x_range(1) ) / nstep;

x(1) = x_range(1);
y(:,1) = y_initial;
yp(:,1) =  feval ( f, x(1), y(:,1) );
%
%  Take one step of the Runge-Kutta-2 method.
%
i = 1;
  xhalf = x(i) + 0.5 * dx;
  yhalf = y(:,i) + 0.5 * dx * yp(:,1);
  yphalf = feval ( f, xhalf, yhalf );

  x(i+1) = x(i) + dx;
  y(:,i+1) = y(:,i) + dx * yphalf;
  yp(:,i+1) = feval ( f, x(i+1), y(:,i+1) );
%
%  After the first step, we proceed with Adams-Bashforth-2.
%
for i = 2 : nstep
  x(i+1) = x(i) + dx;
  y(:,i+1) = y(:,i) + dx * ( 3.0 * yp(:,i) - yp(:,i-1) ) / 2.0;
  yp(:,i+1) = feval ( f, x(i+1), y(:,i+1) );
end