function x = nikola_petrov_gauss_elim ( A, b )

%GAUSS_ELIM   solves the linear system Ax = b using Gaussian 
%             elimination with back substitution
%
%     calling sequences:
%             x = gauss_elim ( A, b )
%             gauss_elim ( A, b )
%     inputs:
%             A   coefficient matrix for linear system 
%                 (A must be square)
%             b   right-hand side vector
%     output:
%             x   solution vector (i.e., vector for which Ax = b)
%
%     NOTE:   this is a demonstration routine - no pivoting strategy 
%             is implemented to reduce the effects of roundoff error

format long e;
[nrow ncol] = size ( A );
if ( nrow ~= ncol )
   disp ('Error: Square coefficient matrix required');
   return;
end;
nb = size ( b );
if ( nrow ~= nb(1) | nb(2) ~= 1 )
   disp ('Error: Size of b incompatible with size of A')
   return;
end;

%
%    Gaussian elimination
%

for i = 1 : nrow - 1
    if ( A(i,i) == 0 )
	   t =  min ( find ( A(i+1:nrow,i) ) + i );
	   if ( isempty(t) ) 
	      disp ('Error: Singular coefficient matrix');
		  return
	   end;
	   temp = A(i,:);     tb = b(i);
	   A(i,:) = A(t,:);   b(i) = b(t);
	   A(t,:) = temp;     b(t) = tb;
	end;
    for j = i+1 : nrow
	    m = -A(j,i) / A(i,i);
		A(j,i) = 0;
		A(j, i+1:nrow) = A(j, i+1:nrow) + m * A(i, i+1:nrow);
		b(j) = b(j) + m * b(i);
    end;
    A
end;
if ( A(nrow,nrow) == 0 )
    disp ('Error: Singular coefficient matrix');
    return
end

%
%    Back substitution
%

x = zeros(nrow,1);
x(nrow) = b(nrow) / A(nrow, nrow);
for i = nrow - 1 : -1 : 1
    x(i) = ( b(i) - sum ( A(i, i+1:nrow) .* x(i+1:nrow)' ) ) / A(i,i);
end;