clear;
a    = 0.;
b    = 1.;
nx   = 55;
h    = 1./(nx+1);

// maillage en x qui contient nx points
x    = [a+h:h:b-h]'

// on construit la fonction phi (donnant la non linéarité)
function y = phi(x)
  y = x.*abs(x);
endfunction

// on construit la fonction phi' (donnant la derivee de la non linéarité)
function y = phip(x)
    //y=1.;
    if x > 0, y = 2*x;
    else, y = -2*x; end
    
endfunction

// on construit le fonction source f(x)
function y = f(x)
    if ((x>0.5)&(x<0.8)) then
        y =10.;
    else y = 0.;
    end
endfunction

// on construit la fonction Fnewton(U,nx) = 0

function Z = Fnewton(U,x,nx,h)
   Z     = zeros(nx,1);  
   Z(1)  = -( phi((U(2)-U(1))/h)  - phi((U(1) -0.00000)/h) )/h - f(x(1));
   Z(nx) = -( phi((0.00-U(nx))/h) - phi((U(nx)-U(nx-1))/h) )/h - f(x(nx));
   for i=2:nx-1
     Z(i) = -( phi((U(i+1)-U(i))/h) - phi((U(i)-U(i-1))/h) )/h - f(x(i));
   end
endfunction
   
// on construit la fonction DFnewton(U,h)
function A = DFnewton(U,x,nx,h)
   h2      = h*h;
   A       = zeros(nx,nx);  
   A(1,1)  = +( phip((U(2)-U(1))/h) + phip((U(1)-0)/h) )/h2;
   A(1,2)  = -( phip((U(2)-U(1))/h))/h2;
   
   A(nx,nx)   = +( phip((0-U(nx))/h) + phip((U(nx)-U(nx-1))/h) )/h2;
   A(nx,nx-1) = -( phip((U(nx)-U(nx-1))/h) )/h2;

   for i=2:nx-1
     A(i,i)   = +( phip((U(i+1)-U(i))/h) + phip((U(i)-U(i-1))/h) )/h2;
     A(i,i+1) = -( phip((U(i+1)-U(i))/h) )/h2;
     A(i,i-1) = -( phip((U(i)-U(i-1))/h) )/h2;
   end
endfunction

// on applique l'algorithme de Newton
err = 1.;
K0  = 50;
U0  = (x-a).*(b-x);
k = 0;

while ((err > 1.E-6)&(k < K0)) do
 FF = Fnewton(U0,x,nx,h);
 AA = DFnewton(U0,x,nx,h);
 YY = AA\FF;
 U1 = U0 - YY;
 err = sqrt(h)*norm(U1-U0)
 U0=U1;
 k   = k+1
 plot(x,U1)
end




