function [L,U,x]=LUDoolittle(A,b)
n=length(A);
L=eye(n);
for k=1:n
    if k==1
        U(k,:)=A(k,:);
        L(k+1:n,k)=A(k+1:n,k)/U(k,k);
    else
        U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
        if k~=n
            L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
        end
    end
end
for k=1:n
    if k==1
        y(k,1)=b(k,1);
    else
        y(k,1)=b(k,1)-L(k,1:k-1)*y(1:k-1,1);
    end
end
for k=n:-1:1
    if k==n
        x(k,1)=y(k,1)/U(k,k);
    else
        x(k,1)=(y(k,1)-U(k,k+1:n)*x(k+1:n,1))/U(k,k);
    end
end
end