function [Ab,x] = Naive_Gaussian(A, b)
% Implements the Gaussian algorithm for solving systems of linear equations
% A coefficient matrix (must be squared matrix)
% b vector of constants
% solve A x = b

[m,n] = size(A);
if  m~= n  % not equal
    error('Matrix A must be square')
end
if length(b) ~= n 
    error('Matrix A and vector b must have same number of rows')
end

Ab = [A b]; % extended coefficient matrix

% Forward elimination (all elements below diagonal set to zero)

for i = 1:n-1  % treating the columns
    for j = i+1:n  % treating each row in i-th column
        p = Ab(j,i)./Ab(i,i); % prefactor
        Ab(j, :) = Ab(j, :) - p.* Ab(i, :)  % Addition of lines
    end 

end

x = zeros(n,1);  % initialization of solution vector
x(n) = Ab(n,n+1)./Ab(n,n);
for i = n :-1:1
    x(i) = (Ab(i, n+1)-Ab(i, i+1:n) *   x(i+1:n)         )./Ab(i,i);
    %                  Block of A  MM  block of already found solutions
    % Note: Sum is included in matrix multiplication! 
end