%% 
% The system is again
% defined by its residual and Jacobian, but this time we implement them as
% a single function.
clear
clc
    
%%
% Our initial guess is the origin. The output has one column per iteration.
x1 = [0;0;0];   % column vector!
x = newtonsys(@nlsystem, x1, 20);
[~, num_iter] = size(x)

%%
% The last column contains the final Newton estimate. We'll compute the
% residual there in order to check the quality of the result.
r = x(:, end)
err = norm(nlsystem(r))

%%title
% Looking at the convergence of the first component, we find a
% subquadratic convergence rate, just as with the secant method.
semilogy(abs(x(1, 1:end-1) - r(1)), 'o-')
grid on
title('Convergence of Newton iterations')
xlabel('n')
ylabel('|x_n(1) - r(1)|')