function y = ABM4PredCorr(f,x,y0)
%
% ABM4PredCorr uses the fourth-order Adams-Bashforth-Moulton predictor-
% corrector formula to solve a first-order ODE in the form y' = f(x,y).
%
%   y = ABM4PredCorr(f,x,y0) where
%
%       f is an inline function representing f(x,y),
%       x is a vector representing the mesh points,
%       y0 is a scalar representing the initial value of y,
%  
%       y is the vector of solution estimates at the mesh points.
%       
py = zeros(4,1);   % Pre-allocate
y(1:4) = RK4(f,x(1:4),y0);   % Find the first 4 elements by RK4
h = x(2) - x(1);

% Start ABM4
for n = 4:length(x)-1,
    py(n+1) = y(n) + (h/24)*(55*f(x(n),y(n))-59*f(x(n-1),y(n-1))+37*f(x(n-2),y(n-2))-9*f(x(n-3),y(n-3)));
    y(n+1) = y(n) + (h/24)*(9*f(x(n+1),py(n+1))+19*f(x(n),y(n))-5*f(x(n-1),y(n-1))+f(x(n-2),y(n-2)));
end