%龍格--庫塔公式
function R=adams(f,a,b,ya,N)
%y'=f(x,y)
%a,b左右端點
%N為迭代步長
%h為步長
%ya為初值
h=(b-a)/N;
T=zeros(1,N+1);
Y=zeros(1,N+1);
T=a:h:b;
Y(1)=ya;
for j=1:N
    k1=h*feval(f,T(j),Y(j));
    k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
    k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
    k4=h*feval(f,T(j)+h,Y(j)+k3);
    Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end
f0=Y(1);
f1=-2*T(2)*Y(2);
f2=-2*T(3)*Y(3);
f3=-2*T(4)*Y(4);
k=4;
    Y1(k+1)=Y(k)+h/24*(-9*f0+37*f1-59*f2+55*f3)
    Y(k+1)=Y(k)+h/24*(f1-5*f2+19*f3+9*feval(f,T(k+1),Y1(k+1)))


R=[T' Y'];