% [x,diffs] = gauss_seidel(A,b,x_0,tol) takes a square matrix A and a vector
% b, a first guess x_0 and a tolerance tol and returns an approximation
% of the solution of Ax=b using Jacobi iteration.
%
% diffs(m) is || x^(m)-x^(m-1)||

function [x,diffs] = gauss_seidel(A,b,x_0,tol)

n = length(b);

% set x to be the first guess:
x = x_0;
% need to define x_old as a temporary copy of x 
x_old = x_0;

ii = 1;
% set test to be large enough that we enter the while loop.
test = 2*tol;
diffs(ii) = test;

while test > tol
  for i=1:n
    x(i) = b(i);
    for j=1:n
% for j less than i, we already know the new values of x.  Use them 
% in defining the new value of x(i).
      if j < i
        x(i) = x(i) - A(i,j)*x(j);       
      elseif j > i
% for j larger than i, we only know the old values of x.  Use them 
% in defining the new value of x(i).
        x(i) = x(i) - A(i,j)*x_old(j);
      end
    end
    x(i) = x(i)/A(i,i);
  end
  test = max(abs(x_old-x));	
  ii = ii + 1;
  diffs(ii) = test;
  x_old = x;
end