function [out,Tj,i] = jacobi(L,U,D,b,x0,n)
x = x0;
Tj = D\(L+U);
cj = D\b;
x(:,2) = Tj*x(:,1)+cj;
i = 2;
while i<=n & ...
        norm(x(:,i)-x(:,i-1),inf)*(1+1/norm(x(:,i),inf))>1e-2;
    i = i+1;
    x(:,i) = Tj*x(:,i-1) + cj;
end
out = x;
return;