function [L,U] = lunaive(A)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LU Factorization (first  %
% principles)              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[rows,cols] = size(A);

if rows ~= cols
   disp('Matrix is not square!')
else
   n = rows;
   L = eye(n);
   U = A;
   
   for p = 1:n-1
      L(p+1:n,p) = U(p+1:n,p)/U(p,p);
      U(p+1:n,p+1:n) = U(p+1:n,p+1:n) - L(p+1:n,p)*U(p,p+1:n);
   end
   U = triu(U);
end
