clc
A=[1 1 3 1;-2 3 5 1;3 0 4 -2;5 -4 2 6];
b=[3 0 -5 4]';
%LU分解法（先進行LU分解）
[l,u]=lu(A);
x0=u\(l\b)
%LU分解法
x1=A\b

n=length(A);
U(1,:)=A(1,:);% 第1行
L(:,1)=A(:,1)/U(1,1);% 第1列
for i=1:n
    L(i,i)=1;% L的對角線為1
end
for k=2:n
    for j=k:n
        U(k,j)=A(k,j);% 第k行
        for r=1:k-1
            U(k,j)=U(k,j)-L(k,r)*U(r,j);
        end
    end
    if k<n
        for i=k+1:n
            L(i,k)=A(i,k)/U(k,k);% 第k列
            for r=1:k-1
                L(i,k)=L(i,k)-L(i,r)*U(r,k)/U(k,k);
            end
        end
    end
end 
% 計算y(i)
for i=1:n
    y(i)=b(i);
    for j=1:i-1
        y(i)=y(i)-L(i,j)*y(j);
    end
end
% 計算x(i)
for i=n:-1:1
    x(i)=y(i)/U(i,i);
    for j=i+1:n
       x(i)=x(i)- U(i,j)*x(j)/U(i,i);
    end
end
x
y   
U
L