function p=steffensen(g,p0,tol,max_iterations)

% Steffensen's method: See section 2.5.
for i=1:max_iterations
    % Store the previous iteration in p0.
    p1 = g(p0);
    p2 = g(p1);
    p = p0 - (p1-p0)^2 / (p2 - 2*p1 + p0);
    if abs(p-p0)<tol
        return
    end
    p0=p;
end
error('Exceeded maximum number of iterations');