function LU = lu_fact(A)
% LU factorization of A: it returns a single matrix LU containing both L
% and U as [L\U].

n = size(A,1);

for j = 1:n-1
    for i = j+1:n
        if A(i,j) ~= 0
            mult = A(i,j)/A(j,j);
            A(i,j+1:n) = A(i,j+1:n) - mult*A(j,j+1:n);
            A(i,j) = mult;
        end
    end
end

LU = A;

end