function p = newton(f, df, p0, tol)
% Solve f(p) = 0 using Newton's method.

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

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