10.3 Heat Conduction in a Sphere with a Nuclear Heat Source
>  restart;
>  A:=r-> 4*Pi*r^2;  area of conduction in a spherical shell
>  Sn := r->Sn0*(1+b*(r/RF)^2);  define the volumetric nuclear energy source: Sn as in eq. 10.3-1 
>  Q:= r->qr(r)*A(r);  eq. 10.3-2 giving the thermal energy flowing across the surface of a sphere of rarius r
>  LHS:=limit((Q(r+dr)-Q(r))/dr,dr=0);  Left hand side of eq. 10.3-5
>  de1:=simplify((LHS=Sn(r)*A(r))/(4*Pi));  Equivalent to eq. 10.3-6

>  de2:=simplify((LHS=0)/(4*Pi)); The DE for the cladding
>  s:=dsolve({de1,qr(0)=finite},qr(r));  solving the ODE and using the BC 10.3-10
>  assign(s);  We must assign the result given by dsolve and
>  qF:=unapply(qr(r),r);  use unapply to make the flux in the fissionable region: qF a function. This agrees with 10.3-12

>  qr:='qr';  We must get rid of the value assigned to qr to use it in the next step.
>  s:=dsolve({de2,qr(RF)=qF(RF)},qr(r));  Now we can find qr in the cladding using 10.3-7 and the BC 10.3-11
>  assign(s);qC:=unapply(qr(r),r);  eq. 10.3-13

>  de3:= -kF*D(TF)(r)=qF(r); de4:=-kC*D(TC)(r)=qC(r);  equations 10.3-14 and 10.3-15: applications of Fourier's Law
>  s:=dsolve({de4,TC(RC)=T0},TC(r)); assign(s);  solve the ODE for the cladding with BC 10.3-19
>  TC:=unapply(TC(r),r);
>  dTC:=r->simplify(TC(r)-T0);  This agrees with eq. 10.3-21 as seen by:
>  dTCbook:=Sn0*RF^2/(3*kC)*(1+3*b/5)*((RF/r)-(RF/RC));
>  simplify(dTC(r)-dTCbook);
>  s:=dsolve({de3,TF(RF)=TC(RF)},TF(r)); assign(s);  solve 10.3-14 with BC 10.3-18
>  TF := unapply(TF(r),r);
>  dTF:=r->simplify(TF(r)-T0);dTFbook:=Sn0*RF^2/(6*kF)*((1-(r/RF)^2)+3*b/10*(1-(r/RF)^4))+Sn0*RF^2/(3*kC)*(1+3*b/5)*(1-RF/RC);differ:=simplify(dTF(r)-dTFbook); The difference TF(r) - T0 we found compared to the result in eq. 10.3-20. 
>  with(plots):
>  tauC:=rho->simplify(dTC(RC*rho)/(Sn0*RF^2/3/kC)); rho is the dimesionless variable: r/RC, tauC is (TC-T0)*3*kc/(Sn0*RF^2)
>  tauC(rho);
>  RF:=alpha*RC;tauC:=unapply(tauC(rho),rho,alpha,b); alpha=RF/RC
>  tauF:=rho->simplify(dTF(RC*rho)/(Sn0*RF^2/3/kC)); tauF is (TF-T0)*3*kc/(Sn0*RF^2)
>  tauF(rho);
>  kF:=beta*kC; beta=kF/kC
>  tauF:=unapply(tauF(rho),rho,alpha,beta,b);
>  p1:=plot(tauF(rho,.75,1.5,.2),rho=0...0.75,title=`(T-T0)*3*kc/(Sn0*RF^2) vs r/RC`): alpha=0.75, beta=1.5, b=0.2, rho goes from 0 to 0.75
>  p2:=plot(tauC(rho,.75,.2),rho=0.75....1.0): alpha=0.75, b=0.2, rho goes from .75 to 1
>  display({p1,p2}); Note the change in slope at rho=0.75
[Maple Plot]
>  ?display
>