function [L,R,fEvals] = Bisection(f,L0,R0,delta,fEvalsMax)
% f is a handle to a continuous function of a single variable.
% L0 and R0 define an interval [L0,R0] and f(L0)f(R0) <= 0.
% delta is the length of an acceptable final interval.
% fEvalsMax is a positive integer >= 2 that indicates the maximum
%    number of f-evaluations allowed.
%
% Bisection is applied to obtain a new bracketing interval [L,R]
% that usually satisfies R-L <= delta. If this is not possible, then
% R-L is as small as possible given the fEvalsMax constraint.

L = L0; fL = f(L); 
R = R0; fR = f(R); 
fEvals = 2;
while (R-L > delta) && (fEvals < fEvalsMax)
    m = (L+R)/2; fm = f(m); 
    fEvals = fEvals+1;
    if fL*fm<=0
        % There is a root in [L,m].
        R  = m; fR = fm;
    else
        % There is a root in [m,R].
        L  = m; fL = fm;
    end
end