function [x, fevals] = bisec(fun, a, b, tol)
%BISEC - calculate the solution of (nonlinear) equation fun(x) = 0
%        using the method of bisection
%        Input: the function handle fun, a and b such that a < b and
%               fun(a)*fun(b) < 0, and a tolerance
%        Output: approximate solution such that abs(xexact - x) <= tol
%                and the number of function evaluations
    fa = fun(a);
    fb = fun(b);
    fevals = 2;
    while (b - a)  > tol
      c = (a + b)/2;
      fc = fun(c);
      fevals = fevals + 1;
      if fc == 0
          break
      end
      if sign(fc)*sign(fa) < 0
          b = c;
          fb = fc;
      else
          a = c;
          fa = fc;
      end
    end
    x = (a + b)/2;
end
