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

D = diag(diag(A));
L = tril(A,-1);
LplusD = L+D;
u = u0;
r0 = b-A*u;
normr0 = norm(r0);
r = r0;
it = 0;

while norm(r)>prec*normr0 && it<itMax
    u = LplusD \ (b-L'*u);
    r = b-A*u;
    it = it+1;
end

