function [root, iterations] = newtons_method(f, f_prime, epsilon, x0)
  ## Find a root of the function `f` with initial guess x0.
  ##
  ## INPUTS:
  ##
  ##   * ``f`` - The function whose root we seek. Must return a column
  ##     vector or scalar.
  ##
  ##   * ``f_prime`` - The derivative or Jacobian of ``f``.
  ##
  ##   * ``epsilon`` - We stop when the value of `f` becomes less
  ##     than epsilon.
  ##
  ##   * ``x0`` - An initial guess. Either 1d or a column vector.
  ##
  ## OUTPUTS:
  ##
  ##   * ``root`` - The root that we found.
  ##
  ##   * ``iterations`` - The number of iterations that we performed
  ##     during the search.
  ##

  if (columns(x0) > 1)
    root = NA;
    return;
  end

  iterations = 0;
  xn = x0;
  f_val = f(x0);

  while (norm(f_val, Inf) > epsilon)
    ## This uses the modified algorithm (2.11.5) in Atkinson.
    ## Should work in 1d, too.    
    delta = f_prime(xn) \ f_val;
    xn = xn - delta;
    f_val = f(xn);
    iterations = iterations + 1;
  end

  root = xn;
end
