adm.m

研一

function Q=abm(f,T,y)
n=length(T);
if n<5
    return;
end;
F=zeros(1,4);
F=feval(f,T(1:4),Y(1:4));
h=T(2)-T(1);

for k=4:n-1
 %预测 校正 
 p=Y(k)+(h/24)*(F*[-9 37 -59 55]');
 T(K+1)=T(1)+h*K;
 F=[F(2) F(3) F(4) feval(f,T(K+1),p)];
 %校正
 Y(K+1)=Y(k)+(h/24)*(F*[1 -5 19 9]');
 F(4)=feval(f,T(K+1),Y(K+1));
end
Q=[T' Y'];