function fnr=demo_newtraph()

    dg = 2;
    h  = 0.01;

    A = zeros(dg,dg);
    b = zeros(dg,1);

    % initial guess
    xo = 1.0; yo = 1.0;
    %xo = 5.0; yo =-5.0;

    % start iterating
    for n=1:20
        fprintf('Step %3d: (x,y)= (%8.4e , %8.4e); (f1,f2)= (%8.4e , %8.4e)\n', n, xo ,yo, f(xo,yo,1), f(xo,yo,2) );
    
        % Jacobian matrix
        A(1,1) = dfdx(xo,yo, 1,h);
        A(1,2) = dfdy(xo,yo, 1,h);
        A(2,1) = dfdx(xo,yo, 2,h);
        A(2,2) = dfdy(xo,yo, 2,h);
        
        % negative of functional value
        b(1) = -f(xo,yo, 1);
        b(2) = -f(xo,yo, 2);

        % solve the linear system
        sol = A\b;
        %[sol,p]=rref([A b]); sol=sol(:,dg+1);
    
        xo = xo+sol(1);
        yo = yo+sol(2);

        %pause;
    end;

    return;

%%%%%%%%%%%%

function f=f(x,y,fid)
    
    if (fid==1) 
        f = x^2 + 4*x - y^2 - 2*y - 1;
    elseif (fid==2)
        f = x^2 + 5*y - 4;
    end
        
    return;

%%%%%%%%%%%%

function dfdx=dfdx(x,y, fid,h)

    dfdx = 1/(2*h) * ( f(x+h,y, fid) - f(x-h,y, fid) );
    return;

%%%%%%%%%%%%

function dfdy=dfdy(x,y, fid,h)

    dfdy = 1/(2*h) * ( f(x,y+h, fid) - f(x,y-h, fid));
    return;

%%%%%%%%%%%%

