function [X, n] = newton(F, JF, X0, tol = 1e-10, maxit = 100) %X = newton(F, JF, X0, tol, maxit) solves the (nonlinear) %system F(X) = 0 using the Newton's iteration with initial %guess X0. (JF is the Jacobi matrix of F.) for n = 1:maxit %Execute one step of Newton's iteration... X = X0 - feval(JF, X0)\feval(F, X0); %... and check if the new approximation is within prescribed tolerance. if(norm(X - X0) < tol) break; end X0 = X; end %A warning in case the last approximation is not within specified tolerance. if(n == maxit) warning("no convergence after maxit iterations") end