学步园

 abm.m

function A = abm(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  - A=[T?? Y??] where T is the vector of abscissas 

%          and Y is the vector of ordinates

n=length(T);

if n>=5

    F=zeros(1,4);

    F=feval(f, T(1:4), Y(1:4));

    h=T(2)-T(1);

    for k=4:n-1

         %Predictor

         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)];

         %Corrector

         Y(k+1)=Y(k)+(h/24)*(F*[1  -5  19  9]');

         F(4)=feval(f, T(k+1), Y(k+1));

    end

    A=[T' Y'];

end

function f=fun1(t,y)

f=t.^2-y;

T=[0,0.05,0.10,0.15,0.20,0.25,0.03];
Y=[1,0.95127058,0.90516258,0.86179202,0,0,0];
f='fun1';
A=abm(f,T,Y)