function [L,U,x]=Crout(A,b)
[n,m]=size(A);
if n==m & det(A)~=0
    for i=1:n;
        L(i,1)=A(i,1);
        U(1,i)=A(1,i)/L(1,1);U(i,i)=1;
    end
    for k=2:n;
        for i=k:n;
            s=0;
            for r=1:k-1;
                s=s+L(i,r)*U(r,k);
            end
            L(i,k)=A(i,k)-s;
        end
        for j=k+1:n;
            s=0;
            for r=1:k-1;
                s=s+L(k,r)*U(r,j);
            end
            U(k,j)=(A(k,j)-s)/L(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