function [root] = bisect(func,Xlow,Xup,XTol)

Flow = feval(func, Xlow);
Fup = feval(func, Xup);

if Flow*Fup > 0
    disp('oops');
    return
end

Xguess = ((Xup + Xlow)/2);
Fguess = feval(func, (Xguess));

while abs((Xup-Xlow)/Xup) > XTol
    if Flow*Fguess > 0
        Xlow = (Xguess);
        Flow = Fguess;
    else
        Xup  = (Xguess);
        Fup = Fguess;
    end
    Xguess = ((Xup + Xlow)/2);
    Fguess = feval(func, Xguess);
end
root = (Xguess);