> restart;

Beginning with equation 11.1-36,

> eq11136:=diff(Thetas(epsilon,tau),tau)=1/epsilon^2*diff(epsilon^2*diff(Thetas(epsilon,tau),epsilon),epsilon);

> with(inttrans):

Now we take the Laplace transform of equation 11.1-36,

> eq11142:=laplace(eq11136,tau,p);

Substituting in the boundary condition,

> eq11142:=subs(Thetas(epsilon,0)=0,eq11142);

We will rename the Laplace transformed variable in order to reduce clutter.

> eq11142:=subs(laplace(Thetas(epsilon,tau),tau,p)=Thetasbar(epsilon),eq11142);

This is now just an ordinary second order differential equation which can be solved in the usual fashion.

> eq11146:=dsolve(eq11142,Thetasbar(epsilon));

Since cos(0)=1, for Thetasbar(0) to be finite (condition 11.1-44), C1 must be equal to 0.

> _C1:=0:eq11146;

This looks fairly similar to equation 11.1-46 The two are nearly equivalent since sin(sqrt(-p))=I*sinh(sqrt(p)). Now work on the fluid half of the problem.

Starting with equation 11.1-40,

> eq11140:=diff(Thetaf(tau),tau)=-3/B*diff(Thetas(epsilon,tau),epsilon);

> eq11145:=laplace(eq11140,tau,p);

Using the boundary condition,

> eq11145:=subs(Thetaf(0)=1,eq11145);

Making some substitutions to make it clearer what we are looking at, and adding that we want to evaluate the derivative at 1,

> eq11145:=subs({laplace(Thetaf(tau),tau,p)=Thetafbar,diff(laplace(Thetas(epsilon,tau),tau,p),epsilon)=D(Thetasbar)(1)},eq11145);

Now using the result earlier (eq11146) to make Thetasbar into a function of epsilon.

> Thetasbar:=epsilon->_C2*sin(sqrt(-p)*epsilon)/epsilon:

And solving for Thetafbar

> eq11147:=Thetafbar=simplify(solve(eq11145,Thetafbar)):eq11147;

Now to satisfy condition 11.1-43

> eq4:=Thetasbar(1)=rhs(eq11147);

We can use this condition to find the unknown constant of integration.

> _C2:=solve(eq4,_C2);

Now look at the equation for Thetafbar again, now that the constant has been determined.

> eq11148:=Thetafbar=simplify(rhs(eq11147));

We can now take the inverse Laplace transformation of both sides of equation 11.1-48.

> eq3:=Thetaf=invlaplace(rhs(eq11148),p,tau);

Unfortunately, this does not seem to be a Laplace transformation that Maple recognizes. I have tried unsuccessfully to get Maple to duplicate the next step in the book, involving a Heaviside partial fraction expansion. Nothing short of doing the expansion by hand and then typing it into Maple seems to work. While there definitely could be some function of Maple that I am not aware of which would solve this problem, as far as I can determine this method is a dead-end in Maple. Perhaps using the partial differential equation solver would be successful, but I imagine that if one wished to recreate Figure 11.1-4 a numeric solution using a method such as finite elements would be the most practical way to proceed.