function xi=NewtonMaehly(alpha,beta,K,x0);
% NEWTONMAEHLY computes zeros of a three term recurrence polynomial 
%   xi=NewtonMaehly(alpha,beta,K,x0) computes by the Newton-Maehly
%   method K zeros of the polynomial defined by the three term
%   recurence with coefficients alpha and beta. The initial guess x0
%   must be bigger than all real zeros. Uses NewtonCorrection.m

xi=[];
for k=1:K,
  y=x0; x=x0*1.1;
  while y<x                                   % Newton double step
    x=y; r=NewtonCorrection(x,alpha,beta);
    s=sum(1./(x-xi(1:k-1)));                  % Maehly correction 
    y=x-2*r/(1-r*s);
  end
  x0=x;
  y=x-r/(1-r*s); x=y+10;                      % single backward step
  while y<x                                   % Newton single steps
    x=y; r=NewtonCorrection(x,alpha,beta);
    y=x-r/(1-r*sum(1./(x-xi(1:k-1))));
  end
  xi(k)=y;  
  k,y                                         % test output 
end
xi=xi(:);
