eqn = D[u[t, r], t] == D[u[t, r], r, r] + D[u[t, r], r]/r;
int = u[0, r] == 1;
bon1 = u[t, 0.1] == 1;
bon2 = u[t, 1] == 1 + t;
sol = NDSolve[{eqn, int, bon1, bon2}, u, {t, 0, 1}, {r, 0.1, 1}];
Plot3D[u[t, r] /. sol, {t, 0, 1}, {r, 0.1, 1}]
註意這裏把0點給挖掉了，因為那裏是奇點，而Mathematica對偏微分數值計算的此類地方要求比較嚴，不去掉的話運算會出錯。這個就是個最簡單的圓形分布的，外邊界的溫度隨時間t線性增加的溫度隨時間變化圖了。非齊次項的引入是完全類似的。讀懂了我的代碼的話，妳就該知道怎麽加了。