function [T,Y] = abm(f,T,Y)
% [T,Y] = abm(f,T,Y)
% Adams-Bashforth-Moulton solution for y' = f(t,y) with y(a) = ya.
% f is the function, input.
% T is the vector of abscissas, input.
% Y is the vector of ordinates, input.
% Remark. The first four coordinates of T and Y
% must have starting values obtained with RK4.m
% T is the vector of abscissas, output.
% Y is the vector of ordinates, output.
n = length(T);
if n<5, break, end;
f0 = feval(f,T(1),Y(1));
f1 = feval(f,T(2),Y(2));
f2 = feval(f,T(3),Y(3));
f3 = feval(f,T(4),Y(4));
h = T(2)-T(1);
h2 = h/24;
a = T(1);
for k = 4:n-1,
p = Y(k) + h2*(-9*f0 + 37*f1 - 59*f2 + 55*f3);
T(k+1) = a + h*k;
f4 = feval(f,T(k+1),p);
Y(k+1) = Y(k) + h2*(f1 - 5*f2 + 19*f3 + 9*f4);
f0 = f1;
f1 = f2;
f2 = f3;
f3 = feval(f,T(k+1),Y(k+1));
end