function [l,u,x,y]=doolittle(a,b)n=length(a); %計算A矩陣的維數for i=1:n %判斷A矩陣的順序主子式是否為非零 w=det(a(1:i,1:i)); if(w==0) error("Matrix must be positive!"); return; endendu(1,1:n)=a(1,1:n); %計算u的第一行元素l(2:n,1)=a(2:n,1)/u(1,1); %計算l的第一列元素for i=1:n %計算l的對角元素 l(i,i)=1;endfor k=2:n for j=k:n sum=0; for s=1:(k-1) sum=sum+l(k,s)*u(s,j); end u(k,j)=a(k,j)-sum; %計算u的其餘元素 end for i=(k+1):n sum=0; for s=1:(k-1) sum=sum+l(i,s)*u(s,k); end l(i,k)=(a(i,k)-sum)/u(k,k); %計算l的其餘元素 endendy(1,1)=b(1); %計算y的第一個元素for i=2:n sumy=0; for j=1:i-1 sumy=sumy+l(i,j)*y(j,1); end y(i,1)=b(i)-sumy; %計算y的其餘元素endx(n,1)=y(n,1)/u(n,n); %計算x的最後一個元素for i=(n-1):(-1):1 sumx=0; for j=i+1:n sumx=sumx+u(i,j)*x(j,1); end x(i,1)=(y(i,1)-sumx)/u(i,i); %計算x的其餘元素end