function x=lusolv(a,b)
%LU分解求解線性方程組aX=b
if length(a)~=length(b)
    error('Error in inputing!')
    return;
end
n=length(a);
[l,u]=lu(a);
y(1)=b(1);
for i=2:n
    z=0;
    for k=1:i-1
        z=z+l(i,k)*y(k);
    end
    y(i)=b(i)-z;
end

x(n)=y(n)/u(n,n);
for i=n-1:-1:1
    z=0;
    for k=i+1:n
        z=z+u(i,k)*x(k);
    end
    x(i)=(y(i)-z)/u(i,i);
end