function    [x,X]=SOR73(A,b,x0,omega,N,TOL);
%function    [x,X]=SOR73(A,b,x0,omega,N,TOL);
%%%%%%%%%%%%%%%%%%% SOR Iterative Algorithm %%%%%%%%%%%%%%%%%%%%%%
%%     
%%     Solve Ax=b by SOR iteration with initial approximation x0 and
%%                parameter omega
%%   INPUTS:
%%     A -  the coefficent matrix of the system
%%     b -  the "Right Hand Side'' of the equation
%%     x0 -  the initial approximation
%%     omega - SOR parameter
%%     N - maximum number of iterations permitted
%%     TOL - tolerance
%%     
%%   OUTPUT:
%%     x - approximate solution or a message to say max number of
%%                                         iterations is exceeded.
%%     X - all the iterations for x (each as a row);

% Step 1
k=1;
n=length(x0);
flag=0;
% Step 2
while  (k<=N & flag==0);
   for i=1:n;
      z=0;
      for j=1:n 
          if  j<i;
           z=z+A(i,j)*x(j);
          end;
          if  j>i;
           z=z+A(i,j)*x0(j);
          end;
      end;
      x(i,1) = (1-omega)*x0(i,1)+(omega*(-z+b(i))/A(i,i));
   end;
   if norm(x-x0)< TOL;
      flag=1;
   end;
   x0=x;k=k+1;
   X=[X;x'];
end;
if flag==0;
   message = 'maximum number of iterations exceeded'
end;
