9.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->Sno*(1+b*(r/Rf)^2);  define the volumetric nuclear 
energy source: Sn as in eq. 9.3-1 
>  Q:= r->q(r)*A(r);  eq. 9.3-2 giving the thermal energy flowing 
across a sphere of rarius r
>  LHS:=limit((Q(r+dr)-Q(r))/dr,dr=0);  eq. 9.3-5
>  de1:=LHS=Sn(r)*A(r);  Equivalent to eq. 9.3-6

>  de2:=LHS=0; The DE for the cladding
>  s:=dsolve({de1,q(0)=finite},q(r));  solving the ODE and using 
the BC 9.3-10
>  assign(s);  We must assign the result given by dsolve and
>  qf:=unapply(q(r),r);  use unapply to make qf a function. 
This agrees with 9.3-12

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

>  de3:= -kf*diff(Tf(r),r)=qf(r); de4:=-kc*diff(Tc(r),r)=qc(r);  equations 9.3-14 and 9.3-15
>  s:=dsolve({de4,Tc(Rc)=T0},Tc(r)); assign(s);  solve 9.3-15 
with BC 9.3-19
>  Tc:=unapply(Tc(r),r);
>  simplify((Tc(r)-T0)/(1+3*b/5));  This agrees with eq. 9.3-21
>  s:=dsolve({de3,Tf(Rf)=Tc(Rf)},Tf(r)); assign(s);  solve 
9.3-14 with BC 9.3-18
>  Tf := unapply(Tf(r),r);
>  dTf:=simplify(Tf(r)-T0);dTfbook:=Sno*Rf^2/(6*kf)*((1-(r/Rf)^2
)+3*b/10*(1-(r/Rf)^4))+Sno*Rf^2/(3*kc)*(1+3*b/5)*(1-Rf/Rc);
differ:=simplify(dTf-dTfbook); The difference Tf(r) - T0 
we found compared to the result in eq. 9.3-20.