Example 12.1-2 Heating of a Finite slab: Part a 
Using Maple for Separation of Variables  to verify eq. 12.1-31
>  restart;
>  eq1:=D[2](T)(y,t)-alpha*D[1](D[1](T))(y,t); Eq. 12.1-3 for one dimensional transient heat conduction in the y direction

>  with(PDEtools,dchange); Used to switch to dimensionless variables.
>  tr:={T=T1-Theta*(T1-T0),y=b*eta,t=b^2*tau/alpha}; Note the first transformation as specified in eq. 12.1-11. This is written so that Theta will be zero ar the two boundaries and Theta will be 1 for t=0.
>  eq2:=dchange(tr,eq1,[Theta(tau,eta),tau,eta]);
>  eq3:=simplify(-eq2*b^2/(alpha*(T1-T0))); eq. 12.1-14

>  pdsolve(eq3); Seeing what the partial differential equation solver does.

That tells us to try separation of variables and gives the ODE that each one must satisfy. _c1 is the separation constant. We will call the function of eta f and the separation constant -c^2 as BS&L does.
>  s:=dsolve({D(D(f))(eta)=-c^2*f(eta),f(1)=0,f(-1)=0},f(eta)); Note that this gives a trivial answer.
>  s:=dsolve({D(D(f))(eta)=-c^2*f(eta),f(-1)=f(1)},f(eta)); This gives a constant times cos(c*eta). Further, if
f(-1)=0, then cos(c) must be zero, since the denominator is -cos(c)^2. All of this is difficult to automate in our solution. Let's go back to
omitting the BC in the solution of our DE.
>  s:=dsolve(D(D(f))(eta)=-c^2*f(eta),f(eta));
>  assign(s);f:=unapply(f(eta),eta);

>  _C1:=0; Since the BC are symmetric in eta, we anticipate f must be also and eliminate the sin term.
>  _EnvAllSolutions:='true'; This will give us all solutions to the next equation when we use solve.
>  eq1:=f(1)=0;eq2:=f(-1)=0;
>  solve(eq1,c);

>  about(_Z1); Maple does not tell us anything specific about _Z1.

_Z1:

  nothing known about this object

>  c:=(n+1/2)*Pi; assume(n,integer); But we know it had better be an integer.

>  eq1;eq2; Then both our BC will be satisfied.

>  f(eta);

>  s:=dsolve(D(g)(tau)=-c^2*g(tau),g(tau)); Moving on to find g(tau)
>  assign(s);g:=unapply(g(tau),tau);

>  g(tau)*f(eta); This agrees with 12.1-25 except _C2 was used as the constant for both f and g.
>  fn:=unapply(f(eta)/_C2,eta,n); The spacial functions without a constant in front.
>  Thetan:=(tau,eta,n)->g(tau)/_C2*fn(eta,n);
>  Thetan(tau,eta,n); Here are the functions we want to sum to get our answer.

>  Theta:=(tau,eta)->sum(D1[n]*Thetan(tau,eta,n),n=0...infinity); We can not call our coefficients D since that implies differentiation.
>  Theta(tau,eta); This agrees with 12.1-26

>  about(n); n is an integer

Originally n, renamed n~:

  is assumed to be: integer

>  IC:=1=Theta(0,eta); Here is the initial condition as shown in 12.1-27.

>  assume(m,integer);
>  int(lhs(IC)*fn(eta,m),eta=-1...1); This is the integral of the LHS*fn of our initial condition.

>  int(rhs(IC)*fn(eta,m),eta=-1...1); It would be nice if Maple would do the same thing to the RHS.

>  int(fn(eta,m)^2,eta=-1...1); We need to take it one step at a time. The fn functions are normalized on the interval.

>  int(fn(eta,n)*fn(eta,m),eta=-1...1); The fn form an orthonormal set of functions on the interval -1 to 1.

>  D1[n]:=int(lhs(IC)*fn(eta,n),eta=-1...1); Thus this gives the coefficients in our series.
>  Theta(tau,eta); This agrees with 12.1-31.

>  tr;

>  Theta(alpha*t/b^2,y/b); This looks even more like eq. 12.1-31 since Theta is (T1-T)/(T1-T0) as seen in the transformation.