clear
clc
%%
% Our initial guess at a root is
% the origin.
x = [0;0;0];   % column vector!

%%
% We need the value (residual) of the nonlinear function, and its Jacobian,
% at this value for x. 
f = nlvalue(x);
J = nljac(x);

%%
% We solve for the Newton step and compute the new estimate. 
s = -(J\f);
x(:,2) = x(:,1) + s;

%%
% Here is the new residual.
format short e
f(:,2) = nlvalue(x(:,2))

%%
% We don't seem to be especially close to a root. Let's iterate a few more
% times. 
for n = 2:7
    f(:,n) = nlvalue(x(:,n));
    s = -( nljac(x(:,n)) \ f(:,n) );
    x(:,n+1) = x(:,n) + s;
end

%%
% We find the infinity norm of the residuals. 
residual_norm = sum(abs(f),1)

%%
% We don't know an exact answer, so we will take the last computed value as
% its surrogate. 
r = x(:,end)