function [L,U] = myLU(A)
% Call : function [L,U] = myLU(A)
% Input : A - nxn matrix
% Output: L - nxn unit lower triangular
%         U - nxn upper triangular


[~,n] = size(A);

L = eye(n);
for i = 1:n-1
    for j = i+1:n
        L(j,i) = A(j,i)/A(i,i);
        A(j,i:n) = A(j,i:n) - L(j,i)*A(i,i:n);
        A(j,i) = 0;
    end
%    keyboard
end

U = A;
return;

     

