solveNewton.m

function x = solveNewton(F, DF, x0, tol)
% uses the Newton method to solve F(x) = 0
% F      function handle, x, F are vectors of the same length
% DF     Jacobian of F (dF/dx)
% x0     start value of iteration
% tol    relative tolerance of result

x = x0;
err = realmax;

while err > tol
  A = DF(x);
  b = -F(x);
  z = A \ b;
  xNeu = x + z;  
  err = norm(xNeu - x)/norm(x);
  x = xNeu;
end