function x = sor(A,b,x0,omega,Nit)

% Solve Ax=b by SOR iteration. 
%   x0: starting value for the iterations. 
%   Nit : Number of iterations. 
%   omega: Relaxation parameter 

n = length(b); 
x = x0;
for k=1:Nit
    for i=1:n
        t = b(i);
        for j=1:i-1
            t = t-A(i,j)*x(j);
        end
        for j=i+1:n
            t = t-A(i,j)*x(j);
        end
        x(i) = (1-omega)*x(i)+omega*t/A(i,i);
    end
end
