function [X, xStarDiff] = GaussSeidel (A, b, x0)
  s = length(b);
  X(1:s, 1) = x0;
  i = 1;
  do
    for j=1:s
      if j==1
        Y(1) = (b(1)-A(1,2:s)*x0(2:s))/A(1,1);
      elseif j == s
        Y(s) = (b(s)-A(s,1:s-1)*(Y(1:s-1))')/A(s,s);
      else
        Y(j) = (b(j)-A(j, 1:j-1)*(Y(1:j-1))' -A(j, j+1:s)*x0(j+1:s))/A(j,j);
      endif
    endfor
    x0 = Y';
    X(1:s, i+1) = x0;
    i = i + 1;
  until (norm(X(:,end-1)-X(:,end),2)<10^-6)

  xStarDiff = norm((A\b)-(X(:,end)),2);

endfunction