% Newton2Ddemo demonstrates Newton's method in two variables
% Finds solution of [x^3 - 3xy^2 = 8,
%                    3x^2y - y^3 = 0].
%                    
x = [-1; 2];     % Initial approximate solution

% Iteration control

xtol = 1e-8; Ftol = 1e-8; maxit = 10;
dx   = 1+xtol; Fnorm = 1+Ftol; m = 0;

% Create header for displayed results

disp('Iter  x(1)                    x(2)                 ||F(x)||    || dx ||')

while (dx > xtol & Fnorm > Ftol & m < maxit)

% Evaluate residual and Jacobian

  r = [x(1)^3-3*x(1)*x(2)^2-8; 3*x(1)^2*x(2)-x(2)^3];
  J = [3*(x(1)^2-x(2)^2)   -6*x(1)*x(2);
       6*x(1)*x(2)          3*(x(1)^2-x(2)^2) ];

% Compute Newton correction 
% (by LU decomposition + fwd & back substitution)
% and apply to x

  s = - (J \ r);
  x = x + s;

% Display progress

  dx = norm(s); 
  Fnorm = norm(r);
  m = m + 1;
  disp(sprintf(' %2d  %21.17f  %21.17f  %8.2e    %8.2e',m,x(1),x(2),Fnorm,dx))
end
disp(sprintf('The computed solution is [%21.17f %21.17f]',x))