function [root,iter,resi,fail] = mueller(p,x0,x1,x2,eps,N)
%
%  Applies Mueller's method to find one root of a polynomial.
%
%  On entry :
%    p          coefficients of a polynomial;
%    x0,x1,x2   three points for the first parabola,
%               as distinct guesses for a root;
%    eps        requested accuracy for the root;
%    N          maximal number of steps.
%
%  On return :
%    root       approximation for the root;
%    iter       number of iterations spent;
%    resi       the residual value equals abs(p(root));
%    fail       if fail == 1, then N steps did not suffice to
%               reach the desired accuracy, otherwise okay.
%
%  Examples :
%    p = rand(1,6);              % p is a random quintic
%    x0 = rand + rand*sqrt(-1);  % and may have complex roots, so we
%    x1 = rand + rand*sqrt(-1);  % choose 3 complex guesses at random
%    x2 = rand + rand*sqrt(-1);
%    [root,iter,resi,fail] = mueller(p,x0,x1,x2,1.0e-12,15);
%    p = poly([1 2 3 4 5]);      % p = (x-1)*(x-2)*(x-3)*(x-4)*(x-5)
%    [root,iter,resi,fail] = mueller(p,x0,x1,x2,1.0e-12,15);
%
fail = 0;
for iter = 1:N
   fx0 = polyval(p,x0);          % fx0 = p(x0)
   fx1 = polyval(p,x1);          % fx1 = p(x1)
   fx2 = polyval(p,x2);          % fx2 = p(x2)
   h0 = x1 - x0; h1 = x2 - x1;   % construct the interpolating parabola
   d0 = (fx1 - fx0)/h0;
   d1 = (fx2 - fx1)/h1;
   a = (d1 - d0)/(h1 + h0);      % a is coefficient of x^2 in parabola
   b = a*h1 + d1;                % b is coefficient of x in parabola 
   c = fx2;                      % c is constant coefficient in parabola
   rad = sqrt(b*b - 4*a*c);      % apply the quadratic formula
   if(abs(b+rad) > abs(b-rad))   % closest root has largest denominator
      den = b + rad;
   else
      den = b - rad;
   end;
   dx = -2*c/den;
   root = x2 + dx;
   fr = polyval(p,root);         % to compute the residual
   resi = abs(fr);
   f = '%3d : %+.15e %+.15e  %.2e  %.2e\n';
   fprintf(f,iter,real(root),imag(root),abs(dx),resi);
   if( (abs(dx) < eps) | (resi < eps) )
      return;                    % successful return
   end;
   x0 = x1; x1 = x2; x2 = root;  % continue one more iteration
end;
fail = 1;                        % failed return