function [t,y] = AB2(t0,y0,t_end,h,fcn)
%
% function [t,y] = AB2(t0,y0,t_end,h,fcn)
%   
% Solve the initial value problem
%    y' = fcn(t,y),  t0 <= t <= b,  y(t0)=y0
% Use the Adams-Bashforth method of order 2 with  
% a stepsize of h. Euler's method is used for the
% value y1. The user must supply a program for the
% right side function defining the differential
% equation. For some name, say deriv, use a first 
% line of the form
%   function ans=deriv(t,y)
% A sample call to AB2 would be
%   [t,z] = AB2(t0,z0,b,delta,'deriv')
%
% Output:
% The routine AB2 will return two vectors, t and y.
% The vector t will contain the node points
%   t(1)=t0, t(j)=t0+(j-1)*h, j=1,2,...,N
% with 
%   t(N) <= t_end,  t(N)+h > t_end
% The vector y will contain the estimates of the
% solution Y at the node points in t.
%
n = fix((t_end-t0)/h) + 1;
t = linspace(t0,t0+(n-1)*h,n)';
y = zeros(n,1);
y(1) = y0;
ft0 = feval(fcn,t0,y0);
y(2) = y0 + h*ft0;      % Euler's method
for i = 3:n
    ft1 = feval(fcn,t(i-1),y(i-1));
    y(i) = y(i-1) + h*(3*ft1-ft0)/2;
    ft0 = ft1;
end

end % AB2