%Adams預報-校正算法
function[x,y]=Adams(fun,x0,y0,N)
h=1/N;
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=x0;
y(1)=y0;
delta(1)=0;
for n=1:N
    x(n+1)=x(n)+h;
    if n<4
        k1=h*feval(fun,x(n),y(n));
        k2=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k1);
        k3=h*feval(fun,x(n)+1/2*h,y(n)+1/2*k2);
        k4=h*feval(fun,x(n)+h,y(n)+k3);
        y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
    else
        k1=h*feval(fun,x(n),y(n));
        k2=h*feval(fun,x(n-1),y(n-1));
        k3=h*feval(fun,x(n-2),y(n-2));
        k4=h*feval(fun,x(n-3),y(n-3));
        yyubao=y(n)+h/24*(55*k1-59*k2+37*k3-9*k4);
        k0=feval(fun,x(n+1),yyubao);
        y(n+1)=y(n)+h/24*(9*k0+19*k1-5*k2+k3);
    end
end
for n=1:N
    delta(n+1)=abs(y(n+1)-x(n+1));
end
plot(x,y,'b')
plot(x,delta,'r')