function [x,y]=Adams_S(ydot_fun, x0, y0, h, N)
% Adams 預報--校正公式，其中
%    ydot_fun ---  一階微分方程的函數
%    x0, y0   ---  初始條件
%    h        ---  區間步長
%    N        ---  區間的個數
%    x        ---  Xn 構成的向量
%    y        ---  Yn 構成的向量

x=zeros(1,N+1); y=zeros(length(y0),N+1); x(1)=x0; y(:,1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    if n<4
        k1=h*feval(ydot_fun, x(n), y(:,n));
        k2=h*feval(ydot_fun, x(n)+1/2*h, y(:,n)+1/2*k1);
        k3=h*feval(ydot_fun, x(n)+1/2*h, y(:,n)+1/2*k2);
        k4=h*feval(ydot_fun, x(n)+h, y(:,n)+k3);
        y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4);
    else
        k1=feval(ydot_fun, x(n), y(:,n));
        k2=feval(ydot_fun, x(n-1), y(:,n-1));
        k3=feval(ydot_fun, x(n-2), y(:,n-2));
        k4=feval(ydot_fun, x(n-3), y(:,n-3));
        ybar=y(:,n)+h/24*(55*k1-59*k2+37*k3-9*k4);
        k0=feval(ydot_fun, x(n+1), ybar);
        y(:,n+1)=y(:,n)+h/24*(9*k0+19*k1-5*k2+k3);
    end
end