function x = Jacobi(A,b,x0)
  
L = -tril(A,-1);
U = -triu(A,1);
D = diag(diag(A));

Nmax = 500;

xe = [-1;1;-1;1];
x = x0;
for j=1:Nmax
  
  x = D\((L+U)*x+b);
  
  if max(abs(x-xe)) < 10^(-4)
    break
  end
    
end
iters = j