function [x,  conv]=mysor(A, b, x, tol, maxit, omega)
% MYSOR - solve Ax=b using SOR 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
             x(j) = (1 - omega)*x(j) + omega*(c(j) - T(j, :) * x); % SOR
        end
        r = A*x - b;
        delta = r'*r;
        conv = [conv delta];
        iter = iter + 1;
    end
end
