% 
% file name: newton_system.m (Newton's method for 2-dimensional problem)
% March 20, 2016  
%

clear all 
clc 
format long

% define functions and their derivatives 
f1   = @(x, y) 4*x^2-y^2;
f2   = @(x, y) 4*x*y^2-x-1;
dxf1 = @(x, y) 8*x;
dyf1 = @(x, y) -2*y;
dxf2 = @(x, y) 4*y^2-1;
dyf2 = @(x, y) 8*x*y;

%
x0 = 0;              % initial guess
y0 = 1;                           
maxiter = 100;      % maximum number of iterations
tolFun  = 1e-12;    % stopping criterion on the function value
tolX    = 1e-12;    % stopping criterion on (x, y) 

%
f1val = f1(x0, y0);
f2val = f2(x0, y0);
[0, x0, y0, f1val, f2val]

% 
for iter = 1:maxiter
    x = x0;
    y = y0;
    
    % jacobian matrix
    jac = [dxf1(x, y) dyf1(x, y);
           dxf2(x, y) dyf2(x, y)];
    
    % solve linear system Ax = b
    h = jac \ -[f1val; f2val];
    
    % Update solutions
    x0 = x + h(1);
    y0 = y + h(2);
    f1val = f1(x0, y0);
    f2val = f2(x0, y0);
    [iter, x0, y0, f1val, f2val] 
    
    % stopping criterions 
    if (abs(x - x0) < tolX && abs(y - y0) < tolX) || (abs(f2val) < tolFun && abs(f1val) < tolFun)
        break
    end
end 
