function r=mnhf_muller(Fun,x,tol,trace,Nmax)
%MNHF_MULLER finds the root of a function "Fun" using Muller scheme
%
% Fun - name of the external function
% x - vector of length 3, (inital guesses)
% tol - error criterion
% trace - print intermediate results 
% Nmax - maximum number of iterations
%
% Usage  mnhf_muller(@poly1,[-0.5 0.75 2.0],1e-8,1,100)
%        poly1 is the name of the external function.
%        [-0.5 0.75 2.0] are the initial guesses for the root.

% Check inputs.
if nargin < 5, Nmax = 100; end
if nargin < 4, trace = 1;  end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 3)
 error('Please provide three initial guesses')
end

frtil=1.0; % Can select any value with |frtil| > tol.
count=1;   % Initialize counter.
while abs(frtil)>tol
 x  = sort(x);         % Ascending order.
 x2 = x(1); x0 = x(2); x1 = x(3); h2 = x0-x2; h1 = x1-x0; gamma = h2/h1;
 f0 = feval(Fun,x0);  % Function evaluations.
 f1 = feval(Fun,x1);
 f2 = feval(Fun,x2);
 % Compute coefficients a, b and c.
 c = f0;
 a = (gamma*f1-f0*(1.0+gamma)+f2)/(gamma*h1^2.0*(1.0+gamma));
 b = (f1-f0-a*h1^2.0)/h1;
 % Evaluate discriminant; if < 0, set to 0.
 discrim = max(b^2.0-4.0*a*c,0.0);
 % Select the appropriate root of the parabola.
 if b >= 0.0
  rtil = x0-2.0*c/(b+sqrt(discrim));
 else
  rtil = x0-2.0*c/(b-sqrt(discrim));
 end
 frtil = feval(Fun,rtil);
 if trace
  fprintf(1,'%3i %12.5f %12.5f\n', count,rtil,frtil);
 end
 % Decide which estimates to keep and which to discard.
 if abs(rtil-x2) < abs(x1-rtil)
  x=[x0 x2 rtil]; % Discard x1.
 else
  x=[x0 x1 rtil]; % Discard x2.
 end
 
 count = count+1; % Increment counter.

 if count > Nmax
  error('Maximum number of iterations exceeded.\n\n')
 end

end
r = rtil;