restart; sol := dsolve(diff(T(r), r) = -a/(r^2), T(r)); T := unapply(rhs(sol), r); sol := solve({T(r1) = T1, T(r2) = T2}, {a, _C1}); assign(sol); normal(T(r)); Q := 4*k*Pi*a; restart; inicon := T(r1) = T1, T(r2) = T2: deq := diff(T(r), r, r) + diff(T(r) ,r)*2/r = 0; sol := simplify(dsolve({deq,inicon}, T(r))); T := unapply(rhs(sol), r); Q := simplify(-4*Pi*r^2*k*diff(T(r, r1, T1, r2, T2),r)); r1 := 0.1: T1 := 400: r2 := 0.4: k := 12: TTs := [100, 200, 300]: for T2 in TTs do; Tpl[T2] := T(r); Qpl[T2] := evalf(Q); od: plot({seq(Tpl[T2], T2 = TTs)}, r = r1..r2); restart; Subs := u = x/(2*a*sqrt(t)): T(x,t):=V(rhs(Subs)); eq := a^2*diff(T(x,t), x, x) = diff(T(x, t), t); eq := subs(x = solve(Subs, x), eq)*t*4; Sol := dsolve(eq, V(u)); T := unapply(rhs(subs(u=rhs(Subs),Sol)), x, t); assume(a >= 0): assume(x >= 0): BounCon := Ts = limit(T(x,t), x = 0, right); IniCon := Tf = limit(T(x, t), t = 0, right); Sol := solve({IniCon, BounCon}, {_C1, _C2}); assign(Sol); T(x,t); evalb(diff(T(x,t),x,x) - diff(T(x,t), t)/a^2=0); Ts := 300: Tf := 285: a := sqrt(0.003): DVEC := [0, 0.01, 0.05, 0.1, 0.2, 0.3]: for d in DVEC do Td[d] := T(d, t) od: plot({seq(Td[d], d = DVEC)}, t = 0..1800, 298..300); plot3d(T(x, t), x = 0..0.5, t = 0..1800, axes = boxed, orientation = [-50, 50]);