function [x, xs, fs] = bisection_root(f, a, b)

% Finds a root of the function f by narrowing the interval [a, b] using
% bisection.  Always maintains the property that sign(f(a)) ~= sign(f(b))
%
%    f -- a function handle to a single-valued function of one variable
%    a -- a starting lower bound
%    b -- a starting upper bound
%
% It must be true that sign(f(a)) ~= sign(f(b)).

if abs(b) < abs(a)
    [a,b] = deal(b,a);
end
    
fa = f(a);
fb = f(b);

n = 0;
while (1)
    n = n + 1;
    if (n > 100), return; end
    
    % eval at midpoint
    c = a + (b - a) / 2;
    fc = f(c);
    fprintf(1, 'c = %#.16g; f(c) = %#.16g\n', c, fc);

    % move an endpoint
    if sign(fc) == sign(fb)
        b = c;
        fb = fc;
    else
        a = c;
        fa = fc;
    end

    % stash values for plotting
    xs(n) = c;
    fs(n) = fc;

    if fc == 0
        x = c;
        break;
    end

    if abs(b - a) <= eps(a)
        x = a;
        break;
    end
    
end

fprintf('Converged in %d iterations\n', n)