--- /dev/null
+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