%% 
% 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)
back_err = norm(nlsystem(r))