function xx = newtonit(x0,nstep)
% Newton solver for n-dimensional system of equations
%   inputs are
%      x0    : starting vector (column of length n)
%      nstep : number of iterations
%   output is
%      xx    : solution vector
n = length(x0);
fprintf('Newton iteration in %g dimensions ...\n',n)
% iteration loop
for k = 1:nstep
[rr,Jmatrix] = fdef(x0);  % compute vector and derivative matrix
itntable(n,k,x0,rr);      % print current solution
dx = - Jmatrix\rr;        % solve
xx = x0 + dx;             % update
x0 = xx;                  % new solution
end
fprintf('done \n')

function itntable(n,k,xx,rr)
if n==2,  % two-dimensional problem
if k==1; % print header
fprintf(' k     x_k      y_k        |f_k|        |g_k|\n'),
end
fprintf(' %g  %8.5f  %8.5f    %6.3e   %6.3e \n', ...
        k-1, xx(1), xx(2),abs(rr(1)),abs(rr(2))), pause(1)
end
return

