%
% Gauss_Seidel.m - a function that computes using SOR Method
%             an approximation xsol of
%             the solution of Ax=b for given
%             a nxn matrix A and a vector b and a starting
%             estimate x^(0) where xsol=x^(k) and
%             ||x^(k)-x^(k-1)||_inf<tol
% Input - A,b,n,x0,tol,kmax,omega
% Output - xsol
%          xn (||x^(k)-x^(k-1)||)
%          flag (0=convergent, 1=divergent)
%          k
%             
function [xsol,xn,flag,k]=SOR(A,b,n,x0,tol,kmax,omega);
D=diag(diag(A));
L=zeros(n);
U=zeros(n);
for i=2:n
    L(i,1:i-1)=-A(i,1:i-1);
    U(i-1,i:n)=-A(i-1,i:n);
end
xv(:,1)=x0;
k=1;
flagct=0;
while flagct==0,
      xv(:,k+1)=inv(D-omega*L)*(((1-omega)*D+omega*U)*xv(:,k)+b);
      xn(k)=norm(xv(:,k+1)-xv(:,k),inf);
      if xn(k)>tol,
          k=k+1;
          if k>kmax,
              flagct=1;
              flag=1;
              disp('SOR Method does not converge in kmax iterations')
          end
      else
          xsol=xv(:,k+1);
          flagct=1;
          flag=0;
          k=k-1;
          disp('SOR Method converges')
      end
end