function [x, xs, fs] = newton_root(f, fprime, x0)

% Finds a root of the function f by iteratively improving the guess x0
%
%    f -- a function handle to a single-valued function of one variable
%    fprime -- a function handle to the derivative of f
%    x0 -- a starting guess

x = x0;

n = 0;
while (1)
    n = n + 1;
    if (n > 100), return; end

    fx = f(x);
    fpx = fprime(x);

    % (diagnostics)
    fprintf(1, 'x = %#.16g; f(x) = %#.16g\n', x, fx);
    xs(n) = x;
    fs(n) = fx;
        
    xprev = x;
    x = x - fx / fpx;

    if abs(x - xprev) < 2*eps(x)
        break;
    end

end

fprintf('Converged in %d iterations\n', n)