这是我把他的代码改成3维的结果，与我用简单初值解出的结果一样。
我原先的解：
本来把一维换成三维只要改一个字节， m=0 换成 m=2, 但Amsel犯了故儒所指出的错误，他把原点放错了地方，所以如果一改他的程序就要crash。其实这么个简单问题Excel都能解，还用什么有限元，真是大炮打蚊子，还打歪了。附能够正确运行的三维代码：
%jiaozi
function [x_jiaozi,t,t]=jiaozi()
% --------------------------------------------------------------
% parameters
rho_out = 1050;
rho_in = 1050;k_out = 0.4;
k_in = 0.5;Cp_out= 3000;
Cp_in = 4000;thickness_out = 1e-3;
thickness_in = 15e-3;Tw = 95;
T_ini = 25;
% --------------------------------------------------------------
% mesh:
x_out = [0.00005 0.0001 0.0005 0.001 0.05 [0.1:0.1:0.9], 0.95, 0.99, 0.995,0.999,0.9995,1]*thickness_out;
x_in = [0.00005 0.0001 0.0005 0.001 0.05 [0.1:0.1:0.9], 0.95, 0.99, 0.995,0.999,0.9995,1]*thickness_in;
x_jiaozi = [0, x_in, thickness_in + x_out];tmax = 600;
t = linspace(0,tmax,100);% --------------------------------------------------------------
m = 2;
sol = pdepe(m,@jiaozi_pde,@(x)T_ini,@jiaozi_bc,x_jiaozi,t);
T = sol(:,:,1);
surf(x_jiaozi,t,T)%--------------------------------------------------------------
function [c,f,s] = jiaozi_pde(x,t,u,DuDx)
if x < thickness_in
k_jiaozi = k_in;
rho_jiaozi = rho_in;
Cp_jiaozi = Cp_in;
else
k_jiaozi = k_out;
rho_jiaozi = rho_out;
Cp_jiaozi = Cp_out;
end
c = rho_jiaozi * Cp_jiaozi / k_jiaozi ;
f = DuDx;
s = 0;
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = jiaozi_bc(xl,ul,xr,ur,t)
pr = ur-Tw; qr = 0;
pl = 0; ql = 1;
end
end