Section 10.4: Heat Conduction with a Viscous Heat Source
>  restart;
>  de:=diff(ex(x),x)=0; eq. 10.4-2
>  s:=dsolve(de,ex(x));
>  assign(s);ex:=unapply(ex(x),x);
>  _C1:=(rho*v^2/2+rho*Hhat)*vx+taudotvx+qx(x); x component of eq. 9.8-6
>  vx:=0;vy:=0;taudotvx:=tauxx*vx+tauxy*vy+tauxz(x)*vz(x); Using Front Cover
>  tauxz:=x->-mu*D(vz)(x);qx:=x->-k*D(T)(x); Newton's and Fourier's Laws
>  ex(x); Left hand side of eq. 10.4-4
>  vz :=x-> (x/b)*vb;  assuming we have a linear velocity profile
>  _C1; eq. 10.4-5
>  s:=dsolve({%=c1,T(0)=T0},T(x));  solve the equation with BC T(0)=T0 and _C1 replaced with c1
>  assign(s); T:=unapply(T(x),x);   
>  BC2:=T(b)=Tb;
>  c1:=solve(BC2,c1);
>  Theta:=xi->simplify((T(xi*b)-T0)/(Tb-T0)); The LHS of eq. 10.4-9 and xi=x/b
>  Theta(xi); 
>  mu:=Br*k*(Tb-T0)/vb^2; Defining the Brinkman No. and replacing the viscosity with it.
>  Theta:=unapply(Theta(xi),xi,Br);
>  Thet:=simplify(Theta(x/b,Br));Thetbook:=((x/b)+(1/2)*Br*(x/b)*(1-(x/b))); difference:=simplify(Thet-Thetbook);  Here is what we found for the left hand side of eq. 9.4-11 compared to the book's answer.
>  with(plots):
>  p1:=plot(Theta(xi,0.5),xi=0...1,color=green,title=`(T-T0)/(Tb-T0)vs x/b with Br=0.5, 2 and 8`):
>  p2:=plot(Theta(xi,2),xi=0...1,color=blue):
>  p3:=plot(Theta(xi,8),xi=0...1,color=red):
>  display({p1,p2,p3});

>