function [L,U,p] = lu_fac(A)
    n = length(A);
    U = A;
    L = zeros(n);
    p = 1:n;
    for i=1:(n-1)
        if (U(i,i) == 0)
            fprintf('Pivot Necessary\n');
        end
        w = U((i+1):n,i)./U(i,i);
        L((i+1):n,i) = w;
        U((i+1):n,:) = U((i+1):n,:) - w*U(i,:);
    end
    L = L + diag(ones(n,1));
end


