function [u,it] = Jacobi(A,b,u0,prec,itMax)

d = diag(A);
LplusLt = A-diag(d);
u = u0;
r0 = b-A*u;
normr0 = norm(r0);
r = r0;
it = 0;

while norm(r)>prec*normr0 && it<itMax
    u = (b-LplusLt*u)./d;
    r = b-A*u;
    it = it+1;
end

