function r=mnhf_bisect(Fun,x,tol,trace)
%MNHF_BISECT finds root of "Fun" using bisection scheme.
%
% Fun - name of the external function
% x - vector of length 2, (initial guesses)
% tol - tolerance
% trace - print intermediate results
%
% Usage  mnhf_bisect(@poly1,[-0.5 2.0],1e-8,1)
%        poly1 is the name of the external function.
%        [-0.5 2.0] are initial guesses for the root.

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

f = feval(Fun,x); % Fun is assumed to accept a vector

if prod(sign(f)) > 0.0
 error('No sign change - no roots')
end

for i = 1:100
 x3 = 0.5*(x(1)+x(2));    % Update the guess.
 f3 = feval(Fun,x3);      % Function evaluation.
 if trace, fprintf(1,'%3i %12.5f %12.5f\n', i,x3,f3); end
 if abs(f3) < tol % Check for convergenece.
  r = x3; 
  return 
 else % Reset values for x(1), x(2), f(1) and f(2).
  if sign(f3*f(1)) < 0.0
   x(2) = x3;
  else
   x(1) = x3;
  end
  f = feval(Fun,x);  % Fun is assumed to accept a vector.
 end
end