function [L,U,x]=Doolittle(A,b)
[n,m]=size(A);
if n==m & det(A)~=0
    for i=1:n;
    U(1,i)=A(1,i);L(i,i)=1;
    L(i,1)=A(i,1)/U(1,1);
    end
    for k=2:n;
        for j=k:n;
            s=0;
            for r=1:k-1;
                s=s+L(k,r)*U(r,j);
            end
            U(k,j)=A(k,j)-s;
        end
        for i=k+1:n;
            s=0;
            for r=1:k-1;
                s=s+L(i,r)*U(r,k);
            end
            L(i,k)=(A(i,k)-s)/U(k,k);
        end
    end
    y(1)=b(1)/L(1,1);
    for i=2:n;
        y(i)=(b(i)-L(i,1:i-1)*y(1:i-1)')/L(i,i);
    end
    x(n)=y(n)/U(n,n);
    for j=n-1:-1:1;
        x(j)=(y(j)-U(j,j+1:n)*x(j+1:n)')/U(j,j);
    end
else 
       disp('此矩陣不為方陣,不能進行lu分解')
end