function p = fixedpoint(g, p0, tol)
% Solve g(p) = p using fixed-point iteration.

fprintf(' n         p          |p-p0|\n');
fprintf('------------------------------\n');

for n = 1:100
    p = g(p0);
    fprintf('%2d  %12.8f  %12.8f\n', n, p, abs(p - p0));
    if abs(p - p0) < tol
        break
    end
    p0 = p;
end
