function [x,k]=Gauss_Seidel_NEW(A,b)
% Resolución de Ax=b por Gauss-Seidel (forma compacta)
%
 D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1);
 n=size(A,1); x=zeros(n,1); k=0; sm=1;
 while sm>0.001
    y = (D-L)\(U*x + b);
    sm=norm(abs(y-x),inf); k=k+1;
    x = y;
 end
end