function x=gem(A,b)
% gem: solution of the linear system Ax=b by gaussian elimination 
%   without pivoting
% Input: A= square non-singular matrix
%        b= right hand side
%
% Output: x= solution of Ax=b
%
[n,m]=size(A);
p=length(b);
if n~= m
    disp('non-square matrix')
    return
end
if p~=m 
    disp('the r.h.s is not compatible with the matrix')
    return
end

% algorithm

for k=1:n
    if A(k,k)==0 
        disp('warning: the submatrix is singular. Pivoting is needed')
        return
    end
    for i=k+1:n
    mik=A(i,k)/A(k,k);
    for j=k+1:n
        A(i,j)=A(i,j)-mik*A(k,j);
     end
       b(i)=b(i)-mik*b(k);
    end
end

x=back_sub(A,b);


function x=back_sub(U,b)
% back_sub: solution of the upper triangular system Ux=b by
%           backward substitutions
% Input: U= upper triangular matrix (non-singular)
%        b= right hand side
%
% Output: x= solution of Ux=b
%
[n,m]=size(U);
p=length(b);
if n~= m
    disp('non-square matrix')
    return
end
if p~=m 
    disp('the r.h.s is not compatible with the matrix')
    return
end
for i=1:n
    if U(i,i)==0
        disp('the matrix is singular')
        return
    end
end

% algorithm
x=zeros(n,1);
for i=n:-1:1
    x(i)=(b(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);
end
