% A=LU decomposition (without pivot)
% where A is nxn matrix,
% L is a unit lower triangular (diagonal entries are 1's)
% and U is upper triangular

function [L,U]=doolittle(A)

n=length(A);

%set L=I and U=A

L=eye(n);  %identity matrix size n
U=A;       %u_ij = a_ij for any ij


for k=1:n-1
    for j=k+1:n
        
        if (U(k,k)==0)       
        disp('cannot factor')
        end

        L(j,k)=U(j,k)/U(k,k);
        U(j,k:n)=U(j,k:n)-L(j,k)*U(k,k:n);

    end %for j loop
end %for k loop




