function[c,iter_hist] = bisection(a,b,fn_tol,error_tol,maxiter,fun)

% This matlab function uses the bisection method to compute
% a solution of the equation f(x)=0 on [a,b].
  iter_hist = [];
  
% Evaluate the functions at the endpoints.
  fa = feval(fun,a);
  fb = feval(fun,b);
  if fa == 0
      disp(' f(a) = 0. ')
      c = a;
      return
  elseif fb == 0
      disp(' f(b) = 0. ')
      c = b;
      return
  end

% Determine whether the bisection method can be used.
  if sign(fa) == sign(fb)
      disp(' sign(f(a)) = sign(f(b)). Choose another interval. ')
      return
  end
  
% The Bisection Method
  for k = 1:maxiter
      e = (b-a)/2;
      c = a + e;
      fc = feval(fun,c);
      fprintf('iter=%d c=%5.5e f(c)=%5.5e e=%5.5e\n',k,c,fc,e)
      iter_hist = [iter_hist; [k, norm(fc), e]];
      if sign(fa) ~= sign(fc) % then f(a)*f(c) < 0
          b = c;
          fb = fc;
      elseif fc == 0
          disp(' f(c) = 0')
          return
      else                    % then f(c)*f(b) < 0
          a = c;
          fa = fc;
      end
      
      % Check the stopping criteria
      if abs(fc) < fn_tol
          disp(' Function value stopping tolerance met.')
          return
      elseif e < error_tol
          disp(' Step size stopping tolerance met.')
          return          
      elseif k == maxiter
          disp(' Max number of iterations reached.')
          return
      end
  end