%LU factorization of square matrix A
%uses naive Gaussian elimination to build L and U

A = [2,4,0;4,5,1;-6,-9,-1];

[n,n] = size(A);      %A is nxn (A must be square)
L = eye(n);           %start with L = nxn identity matrix
U = A;                %start with U = A

%forward elimination on U = A

for k=1:n-1
    for i=k+1:n
        c = U(i,k)/U(k,k);   %c = multiple of row k to subtract from row i
        L(i,k) = c;
        for j=k:n
            U(i,j) = U(i,j)-c*U(k,j);
        end
    end
end

disp('the L and U matrices are...')
L
U
