function [T,Y] = hamming(f,T,Y)
% [T,Y] = hamming(f,T,Y)
% Hamming's 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);
a = T(1);
pold = 0;
cold = 0;
for k = 4:n-1,
pnew = Y(k-3) + 4*h*(2*f1 - f2 + 2*f3)/3;
pmod = pnew + 112*(cold-pold)/121;
T(k+1) = a + h*k;
f4 = feval(f,T(k+1),pmod);
cnew = (9*Y(k) - Y(k-2) + 3*h*(-f2+2*f3+f4))/8;
Y(k+1) = cnew + 9*(pnew-cnew)/121;
pold = pnew;
cold = cnew;
f1 = f2;
f2 = f3;
f3 = feval(f,T(k+1),Y(k+1));
end