MATLAB LU分解法

function [l,u,y,x]=direct_triangle(a,b)
n=max(size(a));
l=eye(n);
u=zeros(n,n);
for k=1:n
    for i=k:n
        u(k,i)=a(k,i);
    end
end
for k=1:n
    for j=1:k-1
        l(k,j)=a(k,j);
    end
end
y=ones(n,1);
x=ones(n,1);
y(1,1)=b(1,1);
for i=2:n
    s=0;
    for k=1:i-1
        s=s+l(i,k)*y(k,1);
    end
    y(i,1)=b(i,1)-s;
end
 x(n,1)=y(n,1)/u(n,n);
 for j=n-1:-1:1
     s1=0;
     for k1=j+1:n
         s1=s1+u(j,k1)*x(k1,1);    
     end   
     x(j,1)=(y(j,1)-s1)/u(j,j);
 end