function x=matrixJ(A,b,n)
D=diag(diag(A));
U=-tril(A,-1);
L=-triu(A,1);
invD=inv(D);
BJ=invD*(L+U);
f=invD*b;
x=zeros(n,1);
rx=BJ*x+f;
while norm(rx-x)>=1.0e-6
    x=rx;
    rx=BJ*x+f;
end
x=rx;
% D=diag(diag(A));  U=-triu(A,1); L=-tril(A,-1);
% B=D\(L+U);   f=D\b;x0=zeros(n,1);
% y=B*x0+f;
% n=1;
% while norm(y-x0)>=1.0e-6
%     x0=y;
%     y=B*x0+f;
%     n=n+1;
% end
% x=y;