function adamm
%
% compute boundary locus for s-step Adams-Moulton method
%
disp('Compute boundary locus for s-step Adams-Moulton method')
s = input('Choose number of step s =');
%
theta0 = 0.0;
theta1 = 2.0*pi;
dtheta = (theta1-theta0)/2000.;
theta  = theta0:dtheta:theta1;
z = complex(cos(theta),sin(theta));
%
alpha = zeros(1,(s+1));
gamm0 = zeros(1,(s+1));
beta  = zeros(1,(s+1));
%
alpha(2) = -1.0;
alpha(1) =  1.0;
rhoz = polyval(alpha,z);
%
gamm0 = 1.0;
for m=2:(s+1)
    gamm0(m) = 0.0;
    for n=1:(m-1)
        gamm0(m) = gamm0(m)-gamm0(n)/(m-n+1);
    end
end
%gamm0
%
% calculate "beta"
% Ref: Henrici p.197
for m=0:s
    beta(m+1) = 0.0;
    for n=m:s
%       # compute binomial coefficient
        bicoef = gamma(n+1)/gamma(m+1)/gamma(n-m+1);
        beta(m+1) = beta(m+1)+bicoef*gamm0(n+1);
    end
    beta(m+1) = (-1)^m*beta(m+1);
end
%beta
%
sigmaz = polyval(beta,z);
%
hbar = rhoz./sigmaz;
%
xmin = -3.;
xmax = 1.;
ymin = -2.;
ymax = 2.;
%
%clf
plot(real(hbar),imag(hbar),...
     [xmin xmax],[0 0],'r-',[0 0],[ymin ymax],'r-')
%
set(gca,'XLim',[xmin xmax], 'YLim',[ymin ymax])
%
axis image
grid on
