function [x,  conv]=mygaussseidel(A, b, x, tol, maxit)
% MYGAUSSEIDEL - solve Ax=b using Gauss-Seidel iterations
    [n, ~] = size(A);
    T = A;
    d = diag(A);
    for i = 1:n
         T(i,i) = 0;
    end
    c = b ./ d;
    T = T ./ d;
    r = A*x - b;
    delta = r'*r;
    conv = delta;
    delta0 = delta;
    iter = 0;
    while (iter < maxit) && (delta > delta0*tol)
        % Jacobi code: x = c - T*x;
        for j = 1:n
             x(j) = c(j) - T(j, :) * x; % Gauss-Seidel
        end
        r = A*x - b;
        delta = r'*r;
        conv = [conv delta];
        iter = iter + 1;
    end
end
