function p = muller(fun, p0, p1, p2, tol, nmax, verbose)

% muller  Uses fixed Mullers method to 
% find a root of the equation  fun(x) = 0 
% with initial values  p0, p1, p2.  
%
% Input:  fun = the left-hand side of the equation  fun(x) = 0
%         p0, p1, p2 = starting values 
%         tol = x tolerance 
%         nmax = maximum number of iterations
%         verbose = 1 to see all iterations, 0 otherwise
%
% Output: r = a root of  fun(x) = 0  possibly complex 
%             with the desired tolerance

if tol < 3*eps       %  The tolerance is too small
  warning('the tolerance is too small\n');
  return;
end

%  The main part of the program starts
%  here

f0 = feval(fun,p0);
f1 = feval(fun,p1);
f2 = feval(fun,p2);

h1 = p1 - p0;
h2 = p2 - p1;
delta1 = (f1 - f0) / h1 ;
delta2 = (f2 - f1) / h2 ;
d = (delta2 - delta1) / (h2 + h1);

if verbose == 1         %  Printing the headings
  fprintf(['\niters         Re(p)              Im(p)    ']); 
  fprintf(['             Re(f(p))            Im(f(p)) \n\n']);
end

for i = 1:nmax

  if verbose == 1   
    fprintf('%3d  ', i);
    fprintf('%18.14f  %18.14f      ', real(p2), imag(p2));
    fprintf('%18.14f   ', real(feval(fun,p2)) );
    fprintf('%18.14f\n', imag(feval(fun,p2)) );
  end

  b = delta2 + h2 * d ;
  DD = sqrt( b * b - 4 * f2 * d );

  if ( abs(b - DD) < abs(b + DD) )
    E = b + DD;
  else 
    E = b - DD;
  end

  h = -2 * f2 / E ;
  p = p2 + h ;

  if ( abs(h) < tol ) 
    return;
  end;

  p0 = p1;
  p1 = p2;
  p2 = p;
  f0 = feval(fun,p0);
  f1 = feval(fun,p1);
  f2 = feval(fun,p2);

  h1 = p1 - p0;
  h2 = p2 - p1;
  delta1 = (f1 - f0) / h1 ;
  delta2 = (f2 - f1) / h2 ;
  d = (delta2 - delta1) / (h2 + h1);
  
end;

warning('root not found within tol after %d iters\n',i);
