function [sol,a,b,k] = bisect(f,a,b,tol)

fa = feval(f,a); fb = feval(f,b);

k=0;
if fa*fb>0
    fprintf('Two endpoints have same sign\n')
    return
end

tol = max(tol,eps);

while abs(b-a)>tol
    k=k+1;
    m = (a+b)/2;
    fm = feval(f,m);
    if fa*fm<=0, b=m;
    else a = m;
    end
end
sol = (a+b)/2;