function [x,k,err]=Newton(B,b,g,x0,eps,kmax)
//- Methode de Newton pour resoudre min(Bx-b,x-g)=0; 

//- Initialisations:
k=0;
x=x0;
err=eps+1;

Fp=spzeros(B);
Id=speye(B); 
PRINT=0

while( k<kmax & err>eps )
  
  k=k+1 

  xold=x;

  //- Definition de F'(x):
  Fp=B;
  test=(B*x-b) - (x-g); 
  j=find(test>0);  Fp(j,:)=Id(j,:);
  F=min(B*x-b,x-g);

  //- Definition nouvel x
  x=x-Fp\F;
 
  //- Estimateur pour convergence
  err=norm(min(B*x-b,x-g),'inf');

  if PRINT
    printf('k=%5i, err=%12.6f, n(x-xold)=%10.6f\n',k,err,norm(x-xold)); 
    //scanf('%c');
  end

end

endfunction