Doolittle分解（matlab代码）

function [L,U] = myLU(A)
    n=length(A);
    L=eye(n,n);
    U=zeros(n,n);
    U(1,1:n)=A(1,1:n);
    L(2:n,1)=A(2:n,1)./U(1,1);

    for r=2:n
        i=r:n;
        U(r,i)=A(r,i)-L(r,1:r-1)*U(1:r-1,i);
        i=r+1:n;
        L(i,r)=(A(i,r)-L(i,1:r-1)*U(1:r-1,r))/U(r,r);
    end
end