function M=milne(f,T,Y) %Input - f is the function entered as a string 'f' % - T is the vector of abscissas % - Y is the vector of ordinates %Remark. The first four coordinates of T and Y must % have starting values obtained with RK4 %Output - M=[T' Y'] where T is the vector of abscissas and % Y is the vector of ordinates % NUMERICAL METHODS: MATLAB Programs %(c) 1999 by John H. Mathews and Kurtis D. Fink %To accompany the textbook: %NUMERICAL METHODS Using MATLAB, %by John H. Mathews and Kurtis D. Fink %ISBN 0-13-270042-5, (c) 1999 %PRENTICE HALL, INC. %Upper Saddle River, NJ 07458 n=length(T); if n F=zeros(1,4); F=feval(f,T(1:4),Y(1:4)); h=T(2)-T(1); pold=0; yold=0; for k=4:n-1 %Predictor pnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]'); %Modifier pmod=pnew+28*(yold-pold)/29; T(k+1)=T(1)+h*k; F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)]; %Corrector Y(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]'); pold=pnew; yold=Y(k+1); F(4)=feval(f,T(k+1),Y(k+1)); end M=[T' Y'];
