function [L,U] = lufactor(A)
% Function to demonstrate the LU factorization of a matrix A
%  such that LU = A
%
% A is nxn matrix
% L is nxn lower triangular matrix
% U is nxn upper triangular matrix
% Syntax: [L,U] = lufactor(A)

[rA,cA] = size(A);  % Determine the size of A

% Input checks
if rA ~= cA
    disp('Error: A is not a square matrix')
    return
else
    n = rA;
end

% Initialize L and U matrices
L = eye(n);
U = A;

disp(' ')
disp('Original Matrix [A]')
disp('-----------------------------')
disp(A)
disp(' ')

disp('L:')
disp(L)
disp('U:')
disp(U)
disp(' ')

for k = 1:n-1
    disp(['Pivot Row ',num2str(k)])
    disp('-----------------------------')
        
    for i = k+1:n
        m = U(i,k)/U(k,k);
        for j = k:n
            U(i,j) = U(i,j)-m*U(k,j);         
        end
        L(i,k) = m;
    end
    disp('L:')
    disp(L)
    disp('U:')
    disp(U)
    disp(' ')   
end