function p = newton(f, df, p0, tol)
    while 1
        p = p0 - f(p0)/df(p0);
        if abs(p-p0) < tol, break; end
        p0 = p;
    end
end
