//  on calcule la solution par une methode de Newton
//
function U = cal_sol_newton(cal_Melm,gradP,U)
//
    eps = .0001; h = .01; imax = 50;

   du = cal_sol_lin(cal_Melm,gradP,U);

   i = 0;
   while norm(du) >= eps & i <= imax
         i = i + 1;
         du =   cal_sol_lin(cal_Melm,gradP,U);
         U = U + du;

        trace_u_1D(U,2);        
        fprintf(%io(2),'i = %3i ; ||du|| = %f',i,norm(du));
   end;

endfunction
