function x = GEdemo(A,b,prec,ptol)
% GEshow  Show steps in Gauss elimination and back substitution
%         No pivoting is used.
%
% Synopsis:  x = GEshow(A,b)
%            x = GEshow(A,b,ptol)
%
% Input:     A,b  = coefficient matrix and right hand side vector
%            prec = number of sig dig precision arithmetic used
%            ptol = (optional) tolerance for detection of zero pivot
%                   Default:  ptol = 50*eps
%
% Output:    x = solution vector, if solution exists

if nargin<4, ptol = 50*eps;  end
[m,n] = size(A);
if m~=n,  error('A matrix needs to be square');  end
nb = n+1;   Ab = [A b];    %  Augmented system
fprintf('\nBegin forward elmination with Augmented system:\n');  disp(Ab);

% --- Elimination
for i = 1:n-1               %  loop over pivot rows
  pivot = chop(Ab(i,i),prec);
  if abs(pivot)<ptol, error('zero pivot encountered');  end
  for k = i+1:n             %  k = index of next row to be eliminated
     Ab(k,i:nb) = chop(Ab(k,i:nb),prec) - chop(chop((chop(Ab(k,i),prec)/pivot),prec)*chop(Ab(i,i:nb),prec),prec); 
     Ab(k,1:nb) = chop(Ab(k,1:nb),prec);  
  end
  fprintf('\nAfter elimination in column %d with pivot = %f\n',i,pivot);
  disp(Ab);
end

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