function [L,U]=LU(A)
n=size(A,1);
L=eye(n);
U=A(1,:);
for k=1:n,
    U(k,k)=A(k,k)-L(k,1:(k-1))*U(1:(k-1),k);
    for j=(k+1):n,
        U(k,j)=(A(k,j)-L(k,1:(k-1))*U(1:(k-1),j));
    end
    for i=(k+1):n
        L(i,k)=(A(i,k)-L(i,1:(k-1))*U(1:(k-1),k))/U(k,k);
    end
end