clear

function s = f(x)
    s(1,1) = -5.*x(1) + 2.*sin(x(1)) + 2.*cos(x(2));
    s(2,1) = 2.*cos(x(1)) + 2.*sin(x(2)) -5.*x(2); 
endfunction

function A = Df(x)
    A(1,1) = -5. + 2.*cos(x(1));
    A(1,2) = -2.*sin(x(2));    
    A(2,1) = -2.*sin(x(1)); 
    A(2,2) = 2.*cos(x(2)) -5.;     
endfunction


function s = g(x)
    s(1,1) = 2/5*(sin(x(1)) + cos(x(2)));
    s(2,1) = 2/5*(cos(x(1)) + sin(x(2))); 
endfunction


xinit = zeros(2,1);
//xinit = [10000;-10000];


// Resolution par Newton 

x = xinit;

r = f(x);
nr = norm(r);
nr0 = nr;

eps = 1.0E-14;
kmax = 100;

k = 1;
res(k) = log(nr/nr0)/log(10);
while (nr/nr0 > eps)&(k<kmax)



    A = Df(x);
    dx = - A\r;
    x = x + dx;
    r = f(x);
    nr = norm(r);

    k = k + 1;    
    res(k) = log((nr+1.0E-20)/nr0)/log(10);       
   
 
end
scf(1)
xtitle('residu')
plot(res,'-xr')

print(%io(2),x,r)


// Resolution par point fixe contractant 

x = xinit;
r = x-g(x);
nr = norm(r);
nr0 = nr;

eps = 1.0E-14;
kmax = 100;


k = 1;
res2(k) = log(nr/nr0)/log(10);
while (nr/nr0 > eps)&(k<kmax)
    
    x = g(x);
    r = x-g(x);
    nr = norm(r);
    
    k = k + 1;
    res2(k) = log(nr/nr0)/log(10);     
    
end

scf(1)
plot(res2,'-xb')

print(%io(2),x,r)


