function [ x, y ] = ab2 ( f, xRange, yInitial, numSteps )
% [ x, y ] = ab2 ( f, xRange, yInitial, numSteps ) uses
% Adams-Bashforth second-order method to solve a system 
% of first-order ODEs dy/dx=f(x,y).
% f = name of an m-file with signature
%    dydx = f(x,y)
% to compute the right side of the ODE as a column vector
%
% xRange = [x1,x2] where the solution is sought on x1<=x<=x2
% yInitial = column vector of initial values for y at x1
% numSteps = number of equally-sized steps to take from x1 to x2
% x = row vector of values of x
% y = matrix whose k-th row is the approximate solution at x(k).

% M. M. Sussman
% $Date: 2013/01/20 23:58:52 $

x=zeros(1,numSteps+1);
x(1) = xRange(1);
h = ( xRange(2) - xRange(1) ) / numSteps;
y(:,1) = yInitial;

k = 1;
  dydx =  f( x(k), y(:,k) );
  xhalf = x(k) + 0.5 * h;
  yhalf = y(:,k) + 0.5 * h * dydx;
  dydxhalf = f( xhalf, yhalf );

  x(1,k+1) = x(1,k) + h;
  y(:,k+1) = y(:,k) + h * dydxhalf;

for k = 2 : numSteps
  dydxold=dydx;
  dydx = f( x(k), y(:,k) );
  x(1,k+1) = x(1,k) + h;
  y(:,k+1) = y(:,k) + h * ( 3 * dydx - dydxold ) / 2;
end
