function [root,iter,resi,fail] = run_mueller(p,eps,N)
%
%  Applies Mueller's method to find all roots of a polynomial.
%
%  On entry :
%    p          coefficients of a polynomial;
%    eps        requested accuracy for the roots;
%    N          maximal number of steps per root.
%
%  On return :
%    root       approximations for all roots, root(k) approximates
%               the k-th root, where k runs from 1 to length(p)-1;
%    iter       iter(k) is the number of iterations spent on root k;
%    resi       resi(k) is the residual value at the k-th root;
%    fail       if fail(k) == 1, then N steps did not suffice to
%               reach the k-th root at the desired accuracy, 
%               otherwise the k-th roots is accurate enough.
%
%  Examples :
%    p = rand(1,6);           % a random quintic
%    [root,iter,resi,fail] = run_mueller(p,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] = run_mueller(p,1.0e-12,15);
%
x0 = rand + rand*sqrt(-1);
x1 = rand + rand*sqrt(-1);
x2 = rand + rand*sqrt(-1);
deg = length(p) - 1;
root = zeros(deg,1);
iter = zeros(deg,1);
resi = zeros(deg,1);
fail = zeros(deg,1);
q = p;
for k = 1:deg-1
  [root(k,1),iter(k,1),resi(k,1),fail(k,1)] = mueller(q,x0,x1,x2,eps,N);
  d = [1 -root(k,1)];
  q = deconv(q,d);            % deflation step
end;
root(deg,1) = -q(2)/q(1);     % last root is just the last quotient
d = [1 -root(deg,1)];         % divisor to compute the residual
[q,r] = deconv(p,d);          % evaluate by deconvolution
resi(deg,1) = abs(r(2));