function [L,U] = simpleludecomp(A)
%SIMPLELUDECOMP - factor matrux A into a product of upper- 
%                 and lower- triangular matrices

% This version of the code uses vectorization 

    [n, ~] = size(A);
    L = eye(n);

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

    U = A;

end
