function xc = bisect(f,a,b,tol)
if sign(f(a))*sign(f(b)) >= 0
end
fa = f(a);
fb = f(b);
while (b-a)/2>tol
    c = (a+b)/2;
    fc = f(c);
    if fc == 0
        break
    end
    if sign(fc)*sign(fb)<0
        a=c;fa=fc;
    else
        b=c;fb=fc;
    end
end
    xc = (a+b)/2;