gauselim.m

function [x, err] = gauselim(A, b)
% Gaussian Elimination routine
%
% 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;

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

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

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

% FORWARD ELIMINATION

for k = 1:n-1

  [A, err] = gepivot(A, k);   % partial pivoting

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

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

% BACK SUBSTITUTION
if A(n,n) == 0
     err = 1; disp('Matrix is singular'); return
end

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