function [c]=bisect(f,a,b,tol)
% Find the root of the equation f(x)=0
% with tolerance tol using bisection method 

assert(a<b);
fa = f(a);
fb = f(b);
assert(fa*fb<0);

iter = 0;

while b-a > tol,
    iter = iter + 1;
    
    c  = (b+a)/2;
    fc = f(c);
    
    fprintf('Iter: %4d, a: %e, f(a): %e, c: %e, f(c): %e, b: %e, f(b): %e\n',...
            iter,a,fa,c,fc,b,fb);
    
    if fc == 0.0  % in reality this should be fc < tol_fun
        break;
    end
    
    if fa*fc < 0
        b  = c;
        fb = fc;
    else
        a  = c;
        fa = fc;
    end
end
% c is our solution