function x=gemcp6_2(A,b,verbose)
%%   function x=gemcp6_2(A,b,verbose);

%%% Gaussian Elimination with Maximal Column Pivoting
%%     and Backward Substitution  On Square Systems
%% INPUT:
%%   A - coefficient matrix of the system of linear equations
%%   b - Right hand side of the system as a column vector (e.g.  
%%   [1;2;3;4;5])
%%   verbose - Put 1 if you want the result of each main step 
%%   (for Homework)
%%           - Put 0 or leave blank if you only want the solution
%% OUTPUT:
%%   x - solution vector as a column vector
%%   'no solution exists' - if pivoting is not possible
%%   IF verbose==1, the intermediate augmented matrices are given

if nargin <3  % if only two arguments given, then verbose is turned off
  verbose=0;
end;

% Make the augmented matrix;
A=[A,b];

%%Determine the size of the system
n=length(A(:,1));
x=zeros(n,1);
flag=0;

%%STEP 1
        NROW=1:n;

%%STEP 2
for i=1:n-1
%%%%%Step 3--5
  if flag==0;
   M=max(abs(A(NROW(i:n),i)));%%% Max in the column on and 
			      %%% below the diagonal

     if M==0   %%% If max is zero, then no solution
        flag=1;
     else  
     I=find(abs(A(NROW(i:n),i))==M);%%% find first place 
				    %%% where max is attained
          %%%%%%%%(interchange rows if necessary)
        p=I(1)+i-1;
          if p>i;
            NCOPY=NROW(i);
            NROW(i)=NROW(p);
            NROW(p)=NCOPY;
           end;
     end;
%%%%%%Step 6--8 (perform the row operations)
       if flag==0;
         for j=i+1:n
           m=A(NROW(j),i)/A(NROW(i),i);
           A(NROW(j),i:n+1)=A(NROW(j),i:n+1)-m*A(NROW(i),i:n+1);
	   clear m;
         end;
         if verbose==1;
               q=int2str(p);
               col=int2str(i);
               ['column','  ',col,',  ','pivot row','  ',q]
                A(NROW,:)
         end;
       end;
   end;
end;
%%%%STEP 9
if flag==0;
   if A(NROW(n),n)==0
    flag=1;
   end;
end;
%%%Step 10  (begin backward substitution)
if flag==0;
  x(n,1)=A(NROW(n),n+1)/A(NROW(n),n);
%%%%Step 11
  for i=n-1:-1:1;
  x(i,1)=(A(NROW(i),n+1)-sum(x(i+1:n,1)'.*A(NROW(i),i+1:n)))/A(NROW(i),i); 
  end;
else
  x='no unique solution may exist';
end;
