function [x,fail] = mnewton(f,df,x,m,eps,N)
%
%  Applies Newton's method to the function f, with derivative df,
%  starting at the point x and looking for a root of multiplicity m.  
%  It produces a sequences of points and stops when the sequence reaches 
%  the point x where |f(x)| < eps or when two consecutive points in the sequence
%  have distance less than eps apart from each other.
%  Failure is reported (fail = 1) when the accuracy requirement
%  is not satisfied in N steps; otherwise fail = 0 on return.
%
%  Example : 
%  >> f = inline('sin(x)^2*(exp(x)-1)')
%  >> df = inline('sin(2*x)*(exp(x)-1)+exp(x)*sin(x)^2')
%  >> [x,fail] = mnewton(f,df,1.0,3,1e-16,10)
%
%  Cluster example : 
%  >> f = inline('x*(x-1e-10)')
%  >> df = inline('2*x-1e-10')
%  >> [x,fail] = mnewton(f,df,2.0,2,1e-16,10)
fprintf('running the method of Newton for a multiple root...\n');
fx = feval(f,x);
for i = 1:N 
   dx = fx/feval(df,x); 
   x = x - m*dx;
   fx = feval(f,x);
   fprintf('  x = %e, dx = %e, f(x) = %e\n', x, dx, fx);
   if (abs(fx) < eps) | (abs(dx) < eps)
      fail = 0;
      fprintf('succeeded after %d steps\n', i);
      return;
   end
end
fprintf('failed requirements after %d steps\n', N);
fail = 1;