function [L, U] = lufact(A)
% LUFACT   LU factorization (demo only--not stable!).
     n = length(A);
     L = eye(n);
     
     for j = 1:n-1
         for i = j+1:n
             L(i,j) = A(i,j)/A(j,j);
             A(i,:) = A(i,:) - L(i,j)*A(j,:);
         end
     end
     U = triu(A);
end