function x = gauss_seidel(M,b,tol,max_itr)

nb = size(M,1);
nf = size(b,2);
assert(size(M,2) == nb);
assert(size(b,1) == nb);
x = zeros(nb,nf);
for k = 1:nf
    fprintf(2,'resolvendo para a funcao %d\n',k);
    z = zeros(nb,1);
    y = zeros(nb,1);
    err = inf;
    itr = 0;
    fprintf(2,'iniciando gauss seidel tol = %12.10f\n',tol);
    while (err > tol) && (itr < max_itr)
        for i = 1:nb
            soma_1 = 0;
            for j = 1:i-1
                soma_1 = soma_1 + M(i,j)*z(j);
            end
            soma_2 = 0;
            for j = i+1:nb
                soma_2 = soma_2 + M(i,j)*y(j);
            end
            assert(M(i,i) ~= 0)
            z(i) = (b(i,k) - soma_1 - soma_2)/M(i,i);
        end
        err = norm(z-y,inf);
        y = z;
        itr = itr+1;
        fprintf(2,'gauss seidel iteracao %d err = %12.10f\n',itr,err);
    end
    x(:,k) = z;
end