时间显示格式一维瞬态常物性导热方程（左右均为第一类边界条件）

%有借鉴也有自己的理解和修改

clc
clear
%参数输入
length = 3;
TL = 3;
TR = 5;
den = 100;
c = 1000;
k = 10;
s = 10;
dt = 100;
n = 10;
maxstep = 500;

%% 参数定义
ae0 = zeros(n+2,1);
aw0 = zeros(n+2,1);
ap0 = zeros(n+2,1);
ap1 = zeros(n+2,1);
b = zeros(n+2,1);
T1 = zeros(n+2,1);
Tall = zeros(n+2,maxstep+1);

%%
%内部网格
dx=length/n;
ae0(3:n,1) = k/dx;
aw0(3:n,1) = k/dx;
ap0(3:n,1) = den*c*dx/dt-ae0(3:n,1)-aw0(3:n,1);
ap1(3:n,1) = den*c*dx/dt;
b(:,1) = s*dx;

%边界网格
%左边界网格
ae0(2,1) = k/dx;
aw0(2,1) = k/(dx/2);
ap0(2,1) = den*c*dx/dt-ae0(2,1)-aw0(2,1);
ap1(2,1) = den*c*dx/dt;

% 右边界网格
ae0(n+1,1) = k/(dx/2);
aw0(n+1,1) = k/dx;
ap0(n+1,1) = den*c*dx/dt-ae0(n+1,1)-aw0(n+1,1);
ap1(n+1,1) = den*c*dx/dt;

%初始条件
T0(2:n+1,1) = 3;
T0(1,1) = 3;
T0(n+2,1) = 5;

%T1时刻
T1(1,1) = 3;
T1(n+2,1) = 5;
Tall(:,1) = T0;
for t = 1:maxstep
    T1(2:n+1,1) = (ae0(2:n+1,1).*T0(3:n+2,1)+aw0(2:n+1,1).*T0(1:n,1)+ap0(2:n+1,1).*T0(2:n+1,1)+b(2:n+1,1))./ap1(2:n+1,1);
    T0 = T1;
    Tall(:,t+1) = T0;
end
Tall = Tall';
subplot(2,1,1);
plot(Tall(:,2:n+1));%行坐标是1,2,3...，纵坐标是从第2列到第n+1列每一行的元素，共有n条曲线
title('每个控制体的温度随时间变化');
subplot(2,1,2);
p = [0,(1/2)*dx:dx:((2*n-1)/2)*dx,length];%n+2个元素的行向量
plot(p,Tall(end,1:end));%[n+2个元素,n+2个元素]得到一组坐标绘制成曲线
title('最终温度分布');