function [x, err] = gausjord(A, b)
% Gauss - Jordan Routine for solution of single system
%
% USAGE:   [x, err] = gauselim(A, b)
%  where   A = square coefficient matrix, size n X n
%          b = single RHS column vector, size n X 1
%          x = solution column vector, size n X 1
%        err = error code:
%                = 0 = no error
%                = 1 = singular matrix detected
%                = 2 = not a square matrix

err = 0;

%set up space for solution vector
x = zeros(size(b));

% check that input is legal
[n,m] = size(A);
if n ~= m
  disp('Error in gauselim: Matrix must be square.')
  err = 2; return
end

% form augmented matrix by adding b to last column of A
m = n+1; A(:,m) = b;

% ELIMINATION

for k = 1:n

  if k < n
     [A, err] = gepivot(A, k);   % partial pivoting
  end

  if err ~= 0
     disp('Matrix is singular'); return
  end

  for i = 1:n
    if i ~= k
       term = A(i,k)/A(k,k);
       for j = k:m
          A(i,j) = A(i,j) - term*A(k,j);
       end
    end
  end

  % normalize
  temp = A(k,k)
  for j = k:m
    A(k,j) = A(k,j)/temp;
  end

end

x = A(:,m);