% bisect.m: Does bisection method on given function.
% Usage: bisect('f', eps, a, b).  eps is tolerance, [a,b] bracket.
function c = bisect(f, eps, a, b)

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

if fa*fb > 0
    disp('[a,b] doesn''t bracket, moron');
    c = NaN;
    return
end

k=1;
disp('iteration a c b');
while b-a>2*eps %This assures midpoint c is within eps of root.
    
    c = (a+b)/2;
    fc = f(c);
    disp([k,a,c,b]); %Comment this out if you don't want to see progress
    
    if fc*fa < 0 %[a,c] is a bracket
        b = c; fb = fc;
    else %[c,b] is the bracket
        a = c; fa = fc;
    end
    k = k+1;
    
end