function [x,k]=Jacobi_NEW_1(A,b)
% Resuelve Ax=b por Jacobi (forma compacta)
%
 d=diag(A);   r=A-diag(d);
 n=length(b); x=zeros(n,1); k=0; sm=1;
 while sm>0.001
    y = (b-r*x)./d;
    sm=max(abs(x-y))/max(abs(y)); k=k+1;
    x = y;
 end
end