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

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

Id=speye(B); 
PRINT=0

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

  xold=x;

  //- Definition de F=F(x) et Fp=F'(x):
  F=min(B*x-b,x-g);
  Fp=Id; 
  zz=(B*x-b) - (x-g); 
  i=find(zz<=0); Fp(i,:)=B(i,:);

  //- 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