function x = ngauss(A, b)
% ngauss Naive Gaussian elimination
% Syntax: x = ngauss(A,b)
% A is the n by n matrix (not changed by ngauss)
% b is the right hand side vector (not changed by ngauss)
% x is the solution (column vector)

[m,n] = size(A);
if m ~= n
   error('Matrix is not square');
end

% forward elimination

% loop over pivot positions (diagonal elements)
for k = 1 : n-1
   
   % loop down pivot column k beginning below pivot
   for i = k+1 : n
      
      % calculate multiplier for row i
      % and subtract multiple of row k from row i
      mult = A(i,k) / A(k,k);
      for j = 1 : n
         A(i,j) = A(i,j) - mult*A(k,j);
      end
      
      % do the same for the right hand side vector
      b(i) = b(i) - mult*b(k);
   end
end

% Back substitution beginning with last unknown x(n)
% and proceding backward until x(1) is found

x = zeros(n,1);
x(n) = b(n) / A(n,n);
for i = n-1 : -1 : 1;
   % compute unknown x(i) in terms of x(i+1) to x(n)
   s = b(i);
   for j = i+1 : n
      s = s - A(i,j)*x(j);
   end
   x(i) = s / A(i,i);
end

end % ngauss