% Solve the nonlinear systems using Newton iteration method.

clear, clc

% initialization
x = [0.2;0.3;0.4];

ite = 0; % number of iteration
fprintf('initial: x =  (%18.16f,%18.16f,%18.16f)^T \n',x(1),x(2),x(3));
dashString = repmat('-', 1, 90);
fprintf('%s\n', dashString);
fprintf('iter.                x                         log(||err||)         order\n');
fprintf('%s\n', dashString);

% a solution
xtrue = [0.7851969330623552256;0.496611392944656396;0.369922830745872357];

% set the tolerance
tol = 1e-10;

err = 0; % store the sequence of errors
% do the iterarion
while 1
    pause  
    ite = ite + 1;
    dx = -Jacobian(x)\nlfuc(x);
    xnew = x + dx;
    err = [err;norm(xnew - xtrue,2)];
    logerr = log10(err(ite + 1));
    if ite > 2
        order = log(err(ite + 1)/err(ite))/log(err(ite)/err(ite - 1));
    else
        order = 0;
    end
    showIter(xnew,ite,logerr,order);
    if norm(dx) < tol
        break
    else
        x = xnew;
    end
end
fprintf('%s\n', dashString);
fprintf('Matlab: \n')
format long
x = xtrue

function output = nlfuc(x)
% The nonlinear function
output = zeros(3,1);

output(1) = x(1)^2 + x(2)^2 + x(3)^2 - 1;
output(2) = 2*x(1)^2 + x(2)^2 - 4*x(3);
output(3) = 3*x(1)^2 -4*x(2) + x(3)^2;

end

function J = Jacobian(x)
% The Jacobi matrix of the nonlinear function

J = [ 2*x(1), 2*x(2), 2*x(3);
      4*x(1), 2*x(2), -4;
      6*x(1), -4, 2*x(3)];
  
end

function showIter(x,ite,logerr,order)
% Display the information
if ite ~= 0
    for i = 1:3
        if i == 2
            iteString = sprintf('%2d',ite);
            showString = sprintf('%23.16f',x(i));
            if ite > 2
                infoString = [sprintf('%16.8f',logerr) ...
                    sprintf('%12.3f',order)];
            else
                infoString = sprintf('%16.8f',logerr);
            end
            fprintf([iteString showString infoString '\n']);
        else
            fprintf('%26.16f\n',x(i));
        end
    end
else
end
fprintf('\n');
end