% solves f(x) = 0 in [a,b], under the assumption that f(a)*f(b)<0
function [l,r,vl,vr, n] = bisect(a, b, delta, f)
N=10^4;
if a>b
    c=a;
    a=b;
    b=c;
end
l=a;
r=b;
vl=f(l);
vr=f(b);
for i=1:N
    m=(r+l)*0.5;
    vm=f(m);
    if vl*vm<0
        r=m;
        vr=vm;
    else
        l=m;
        vl=vm;
    end
    if (r-l)<delta
        break
    end
end
n=i;