% Gaussian elimination with backward substitution
%This code might not work (on purpuse, haha). It is only meant to be an illustration
% for the students.
% INPUT n; A(i,j) with 1 <= i <= n,  1 <= j <= n+1
% OUTPUT x(1),x(2),...x(n) or error message
% AUTHOR: H. Ceniceros
% LAST MODIFIED: 08/03/09
ich=0; %counter for row exchanges
TRUE  = 1;
FALSE = 0;
OK=TRUE;
ZERO=1.0e-16 ;
% Step 1
  for i=1:n-1; % do steps 2-4
      % Step 2
      ip=i;
      while abs(A(ip,i)) <= ZERO & ip <= n 
         ip = ip + 1;
      end;
      if ip > n 
           OK=FALSE;
           break;
      end
      % Step 3
      if ip = i   %switch rows, yeah right!
         for j = 1:n+1;
             temp = A(i,j);
             A(i,j) = A(ip,j);
             A(ip,j) = temp1;
         end;
         ich=ich + 1  %count row interchanges
      end;
      % Step 4
      for j=i+1:n ;   % Do steps 5 and 6
         % Step 5
         xm = A(j,i)/A(i,i);
         % Step 6
         A(j,i)=0;
         for k=i+1:n+1;
             A(j,k) = A(j,k) - xm * A(i,k);
         end;
     end;
A(1:n,1:n)
 end;
% Step 7
 if abs(A(n,n))  <= ZERO | OK == FALSE 
    disp('no unique solution exists');
    OK=FALSE;
    break
 end;
if OK == TRUE
  % Step 8 (Start backward substitution)
    x(n) = A(n,n+1) / A(n,n);
  % Step 9
    for k = 1:n-1;
        i = n-k;
        sum = 0;
        for j = i+1:n;
            sum = sum +  A(i,j) * x(j);
        end;
         x(i) = (A(i,n+1)- sum) / A(i,i);
     end;
     % Step 10
      disp('procedure completed succesfully');
     x
end;