function x = GE(A,b)
% Basic Gaussian Elimination Routine Adapted from Book's Code
% Note that no pivoting or checking for 0's on the diagonal is done. 
%
% Function call:  x = GE(A,b)
%
% Input:     A,b  = coefficient matrix and right hand side vector
%
% Output:    x = solution vector

[m,n] = size(A);
if m~=n,  error('A matrix needs to be square');  end
nb = n+1;   Ab = [A b];    %  Augmented system


% --- Elimination
for i = 1:n-1               %  loop over pivot rows
  %Choose pivot more intelligently by examining all the entries  in this
  %column on or below the diagonal
  pivot = Ab(i,i);
  for k = i+1:n             %  k = index of next row to be eliminated
     Ab(k,i:nb) = Ab(k,i:nb) - (Ab(k,i)/pivot)*Ab(i,i:nb); 
  end
end

% --- Back substitution
x = zeros(n,1);           %  preallocate memory for and initialize x
x(n) = Ab(n,nb)/Ab(n,n);
for i=n-1:-1:1
  x(i) = (Ab(i,nb) - Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);
end