function [x,res] = weierstrass(p,x,eps,N)
%
%  [x,res] = weierstrass(p,x,eps,N)
%     applies the method of Weierstrass on the polynomial p,
%     starting at the points in the vector x.
%  The iteration stops until the accuracy eps is reached.
%  No more than N steps are performed.
%
%  Required: the polynomial is monic, i.e.: p(1) = 1.
%
%  On entry:
%     p      coefficient vector of polynomial, p(1) = 1;
%     x      starting values for the roots;
%     eps    accuracy requirement on the roots;
%     N      maximum allowed number of iterations.
%
%  On return:
%     x      new vector with approximations for the roots;
%     res    residual, magnitude of p evaluated at x,
%            if (res <= eps) then success, otherwise failure.
%
%  Example :
%     p = rand(1,11); p = p/p(1);
%     x = rand(1,10) + rand(1,10)*i;
%     [x,res] = weierstrass(p,x,1.e-12,100)
%
d = length(x);
fprintf('  real(x)               imag(x)                 dx        f(x) \n');
y = polyval(p,x);
for i = 1:N
   dx = ones(size(x,1),size(x,2));
   for j = 1:d
     for k = 1:d
       if (k ~= j) dx(j) = dx(j)*(x(j) - x(k)); end;
     end;
     dx(j) = y(j)/dx(j);
   end
   x = x - dx;
   y = polyval(p,x);
   fprintf('results of step %d:\n', i);
   for j = 1:d
     fprintf('%22.15e %22.15e ', real(x(j)), imag(x(j)) );
     fprintf(' %.2e %.2e\n', abs(dx(j)), abs(y(j)) );
   end;
   res = norm(y,inf);
   if (abs(dx) < eps) | (res < eps)
      fprintf('succeeded after %d steps\n', i);
      return;
   end
end
fprintf('failed requirements after %d steps\n', N);