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

d = diag(A);

u = u0;
r = b-A*u;
normr = norm(r);
normr0 = normr;
k = 0;

while normr/normr0 > prec
    u = u + r./d;   %  u = u + inv(D)*(b-A*u)
    normr_prev = normr;
    r = b - A*u;
    normr = norm(r);
%     normr/normr_prev, normr/normr0, pause
    k = k+1;
end