function p = bisection(f, a, b, tol)
% Solve f(p) = 0 using the bisection method.

fprintf(' n         a             b             p           f(p)\n');
fprintf('----------------------------------------------------------\n');

for n = 1:100
    p = (a + b) / 2;
    fprintf('%2d  %12.8f  %12.8f  %12.8f  %12.8f\n', ...
        n, a, b, p, f(p));
    if p - a < tol
        break
    end
    if f(a) * f(p) > 0
        a = p;
    else
        b = p;
    end
end
