以向前差分的方法解决二维传热问题的matlab程序

clear; close all; L=100; M=20; timeMax=50000; Currtime = 100; DL=[0:L/M:L]; T = ones(M,M,timeMax)*30; T(1:M,1,1:1:timeMax)=100; T(1:M,M,1:1:timeMax)=100; T(1,1:M,1:1:timeMax)=50; T(M,1:M,1:1:timeMax)=50; f=0.1; for t=2:timeMax for i=2:M-1 for j=2:M-1 T(i,j,t)=f*T(i-1,j,t-1)+f*T(i,j-1,t) +(1-4*f)*T(i,j,t-1) +f*T(i+1,j,t-1)+f*T(i,j+1,t-1); end; end; end; [x,y]=meshgrid(1:M,1:M); mesh(x,y,T(1:M,1:M,Currtime)); title(['2D Temperature Field, t=',num2str(Currtime-1),'s']); xlabel('x length(L/M)'); ylabel('y length(L/M)'); zlabel('Temperature(^oC)');