function [L,U] = Doolittle(A)
% Compute the Doolittle LU factorization of the matrix A
%
% INPUT: A (Matrix to factorize.)

% OUTPUT: L (Unit lower triangular matrix of the factorization.)
%         U (Upper triangular matrix of the factorization.)

% Initialize
n = length(A);
L = eye(size(A));
U = zeros(size(A));

for k = 1:n
    
    % Generate row k of U
    U(k,k:n) = A(k,k:n);
    for j = 1:k-1
        U(k,k:n) = U(k,k:n) - L(k,j)*U(j,k:n);
    end
    
    % Generate column k of L
    L(k+1:n,k) = A(k+1:n,k);
    for j = 1:k-1
        L(k+1:n,k) = L(k+1:n,k) - L(k+1:n,j)*U(j,k);
    end
    L(k+1:n,k) = L(k+1:n,k)/U(k,k);       
end


