function [L,U] = lufact(A)
% LUFACT   LU factorization without pivoting (demo only--not stable!).
% Input:
%   A    square matrix
% Output:
%   L,U  unit lower triangular and upper triangular such that LU=A

    n = length(A);
    L = eye(n);   % ones on diagonal

    % Gaussian elimination
    for j = 1:n-1
        for i = j+1:n
            L(i,j) = A(i,j) / A(j,j);   % row multiplier
            A(i,j:n) = A(i,j:n) - L(i,j)*A(j,j:n);
        end
    end

    U = triu(A);
end