function [L,U]=lufact(A)
%
[n,m]=size(A);

if n~= m
    disp('non-square matrix')
    return
end

% algoritmo

for k=1:n
    if A(k,k)==0 
        disp('warning, the submatrix is singular')
        return
    end
    for i=k+1:n
    A(i,k)=A(i,k)/A(k,k);
    for j=k+1:n
        A(i,j)=A(i,j)-A(i,k)*A(k,j);
     end
    end
end
L=tril(A,-1)+eye(n);
U=triu(A);
