% x = gauss_elim(A,b,i_test) takes a square matrix A and a vector b and 
% returns the solution of Ax=b.  If i_test = 0 then there is pivoting,
% otherwise there is none.  NOTE: the routine writes over A and b, so
% you can't use it twice in a row!

function x = gauss_elim(A,b,i_test)

[n,n] = size(A);

for k=1:n-1
  if i_test == 0
     [A,b] = pivot_98(k,A,b);
  end 
  for i = k+1:n
    m(i,k) = A(i,k)/A(k,k);
    for j=k+1:n
      A(i,j) = A(i,j) - m(i,k)*A(k,j);
    end
    b(i) = b(i) - m(i,k)*b(k);
  end
  A;
end

% now the matrix A is in an upper triangular form.  Ignore the numbers
% below the diagonal, they should be zero, I just didn't force them to
% be.

x(n) = b(n)/A(n,n);
for i=1:n-1
  x(n-i) = b(n-i);
  for j=n-i+1:n
    x(n-i) = x(n-i) - A(n-i,j)*x(j);
  end
  x(n-i) = x(n-i)/A(n-i,n-i);
end

x = x';