非定常の場合について、要は、

T'[k]=F0*[(1+(Δr)/(2*r[k])*T(k+1)+(1/(F0)-2)*T[k]+(1-Δr/(2*r[k]))*T[k-1]+{(Δr)^2}/k*H]
T'[0]=4*F0*{T[1]+(1/(4*F0)-1)*T[0]+{(Δr)^2/(4*k)*H} //境界条件
H=32/(π^2)*10^10 [J/m3・s]

T[k] : K番目における温度
T[k+1] : (k+1)番目における温度
F0 : フーリエ数。Fo=k*Δt/(ρ*Cp*(Δr)^2)
Δr : 分割した時のメッシュの厚さ
r[k] : 円の中心からk番目までの長さ

を解けばいいということになっているのですが、微分方程式を解く所に間違いがありそうです。
添削よろしくお願いしますm(_ _)m