function [T,Y] = milne(f,T,Y)
% [T,Y] = milne(f,T,Y)
% Milne-Simpson 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;
yold = 0;
for k = 4:n-1,
  pnew = Y(k-3) + 4*h*(2*f1 - f2 + 2*f3)/3;
  pmod = pnew + 28*(yold-pold)/29;
  T(k+1) = a + h*k;
  f4 = feval(f,T(k+1),pmod);
  Y(k+1) = Y(k-1) + h*(f2 + 4*f3 + f4)/3;
  pold = pnew;
  yold = Y(k+1);
  f1 = f2;
  f2 = f3;
  f3 = feval(f,T(k+1),Y(k+1));
end
