Heat equation in plane polar coordinates 
  
*
  
*A tall circular cylinder with radius $R$ is in temperature equilibrium with the surroundings which has the temperature $0$.
  
*Heat is produced uniformly from time $t=0$ with the source density $q$.
  
*Determine the temperature distribution over a cylindrical cross section. Disregard from any transition resitance at the surface.
  ( Freely translated exercise )
  

The heat equation together with boundary/initial conditions should look like
$$\cases{\dfrac{\partial u}{\partial t} - a\Delta u = \dfrac{aq}{\lambda} & (1)\\u(R,t) = 0, \: t>0 & (2)
 \\u(r,0) = 0, \:0<r<R & (3)} $$
 where $a$ is the diffusivity constant and $\lambda$ the heat conductivity constant. We assume cylindrical symmetry and that $u$ is bounded. Let's define the operator $\mathcal{A} = - \dfrac{1}{r}\partial_r(r \partial_r)$ with domain $\mathcal{D_A} = \{v \in C^2(0,R) \: | \: v(R) = 0, v(r) \leq C \}$ for some constant $C$. Solving the ODE $\mathcal{A}v = \lambda v$ (not same $\lambda$ as above) together with the boundary condition gives the eigenfunctions $\varphi_k = J_0(\sqrt{\lambda_k}r)$ with corresponding eigenvalues $\lambda_k = (\alpha_{0k}/R)^2$ where $\alpha_{0k}$ are the zeros to the Bessel function $J_0$, $k=1,2,...$ We now assume that $u(r,t) = \sum_{k=1}^{\infty} u_k(t) \varphi_k(r)$, which inserted into $(1)$ yields the equation $$\sum_{k=1}^{\infty} [u_k'(t)+a\dfrac{\alpha_{0k}^2}{R^2}u_k(t)]\varphi_k(r) = \dfrac{aq}{\lambda}$$ A homogenous solution is $\sum_{k=1}^{\infty}  c_k e^{-at \alpha_{0k}^2/R^2}$, but since $u(r,0) = 0 \Rightarrow c_k=0, \forall k$. Now to get the particular solution we run into... 
The problem 
Normally, we (well at least I) would project the right side on the ON-basis $\{\phi_k \}_{k=1}^{\infty}$ but this would involve integrating the stuff like $J_0$ and $J_0^2$ which may or may not be possible. However, the answer sheet provides the answer
\begin{align}
u(r,t) & = \dfrac{qR^2}{4\lambda}(1-(\dfrac{r}{R})^2-8 \sum_{n=1}^{\infty} \dfrac{J_0 (\alpha_{0n}r/R)}{\alpha_{0n}^3J_1(\alpha_{0n})} e^{-\alpha_{0n}^2ta/R^2})
\\[3mm] & = \dfrac{2qaR^2}{\lambda} \sum_{n=1}^{\infty} \dfrac{1}{\alpha_{0n}^3J_1(\alpha_{0n})}(1-e^{-\alpha_{0n}^2at/R^2})J_0(\alpha_{0n}r/R)
\end{align}
which suggests there is a simple approach. Does anyone see it?
 A: I don't see where anything simple is begin applied, just the same old tough-ass Sturm-Liouville song and dance. Let's see, the heat flow vector is
$$\vec q=-k\vec\nabla T$$
and then
$$\vec\nabla\cdot\vec q=-k\nabla^2T=-\rho C_p\frac{\partial T}{\partial t}+\dot q$$
So that's why we started with
$$\nabla^2T-\frac{\rho C_p}k\frac{\partial T}{\partial t}=-\frac{\dot q}k$$
The first step was to take care of that inhomogeneous part, $-\frac{\dot q}k$. Since it's time- and azimuthally-independent, we try to solve
$$\frac1r\frac{d}{dr}\left(r\frac{dT_p}{dr}\right)=-\frac{\dot q}k$$
Integrating,
$$r\frac{dT_p}{dr}=-\frac{\dot q}{2k}r^2+C_1$$
$$T_p=-\frac{\dot q}{4k}r^2+C_1\ln r+C_2$$
Since $T_p$ should be finite as $r\rightarrow0$, $C_1=0$. Then $T_p(R)=-\frac{\dot q}{4k}R^2+C_2=0$, so
$$T_p(r)=\frac{\dot q}{4k}(R^2-r^2)$$
So now we just have to solve the homogeneous equation
$$\nabla^2T-\frac{\rho C_p}k\frac{\partial T}{\partial t}=0$$
And add the particular solution. Given azimuthal and axial independence we can write $T(r,t)=\rho(r)\tau(t)$ and separate variables as
$$\frac1{r\rho}\frac{\partial}{\partial r}\left(r\frac{\partial\rho}{\partial r}\right)=\frac{\rho C_p}{k}\frac{\partial\tau}{\partial t}=-\lambda^2$$
The $\rho$ equation is
$$r^2\frac{d^2\rho}{dr^2}+r\frac{d\rho}{dr}+\lambda^2r^2\rho=0$$
To satisfy the boundary conditions that $\rho(R)=0$ and is also finite at $r=0$ we require $\rho_k(r)=J_0\left(\frac{\lambda_kr}R\right)$ where $J_0(\lambda_k)=0$.  
The $\tau$ equation reads
$$\frac{d\tau_k}{dt}=-\frac{\lambda_kk}{\rho C_p}\tau_k$$
Oh dear, $k$ for both index and thermal conductivity TT. So
$$\tau_k(t)=e^{-\frac{\lambda_kk}{\rho C_p}t}$$
Adding to the particular solution we have
$$T(r,t)=\frac{\dot qR^2}{4k}\left\{1-\frac{r^2}{R^2}+\sum_{k=1}^{\infty}a_kJ_0\left(\frac{\lambda_kr}R\right)e^{-\frac{\lambda_kk}{\rho C_p}t}\right\}$$
Now comes the Sturm-Liouville stuff. I am using Murray R.Spiegel's Schaum's Outline on Fourier Analysis for these formulas I don't carry around in my head. The Bessel functions have normalization
$$\int_0^1xJ_n^2(\lambda_kx)dx=\frac12\left[J_n^{\prime2}(\lambda_k)+\left(1-\frac{n^2}{\lambda_k^2}\right)J_n^2(\lambda_k)\right]$$
Here, $n=0$ and $J_0^{\prime}(x)=-J_1(x)$, $J_0(\lambda_k)=0$ is a boundary condition and using the orthogonality of the Bessel functions
$$\int_0^RrJ_0\left(\frac{\lambda_kr}R\right)J_0\left(\frac{\lambda_jr}R\right)dr=\frac{R^2}2J_1^{2}(\lambda_k)\delta_{jk}$$
So we just have to integrate on both sides like this to pick out the Fourier coefficients $a_k$. We need the identity
$$\frac d{dx}\left[x^nJ_n(x)\right]=x^nJ_{n-1}(x)$$
$$\int_0^RrJ_0\left(\frac{\lambda_kr}R\right)=\frac{R^2}{\lambda_k^2}\int_0^{\lambda_k}xJ_0\left(x\right)dx=\left.\frac{R^2}{\lambda_k^2}xJ_1(x)\right|_0^{\lambda_k}=\frac{R^2}{\lambda_k}J_1(\lambda_k)$$
$$\begin{align}\int_0^Rr\cdot r^2J_0\left(\frac{\lambda_kr}R\right)dr&=
\frac{R^4}{\lambda_k^4}\int_0^{\lambda_k}x^3J_0(x)dx\\
&=\frac{R^4}{\lambda_k^4}\left\{\left.x^2\cdot xJ_1(x)\right|_0^{\lambda_k}-2\int_0^{\lambda_k}x^2J_1(x)dx\right\}\\
&=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-\left.2x^2J_2(x)\right|_0^{\lambda_k}\right\}\\
&=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-2\lambda_k^2J_2(\lambda_k)\right\}\\
&=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-2\lambda_k^2\left(\frac2{\lambda_k}J_1(\lambda_k)-J_0(\lambda_k)\right)\right\}\\
&=R^4\left(\frac1{\lambda_k}-\frac4{\lambda_k^3}\right)J_1(\lambda_k)\end{align}$$
So all this integration amounts to
$$\frac{R^2}{\lambda_k}J_1(\lambda_k)-\frac1{R^2}R^4\left(\frac1{\lambda_k}-\frac4{\lambda_k^3}\right)J_1(\lambda_k)+a_k\frac{R^2}2J_1^{2}(\lambda_k)=0$$
To satisfy the boundary condition that $T(r,0)=0$. So that tells us both that
$$a_k=-\frac8{\lambda_k^3J_1(\lambda_k)}$$
and that
$$1-\frac{r^2}{R^2}=\sum_{k=1}^{\infty}\frac8{\lambda_k^3J_1(\lambda_k)}J_0\left(\frac{\lambda_kr}R\right)$$
So we can put together the solution
$$T(r,t)=\frac{\dot qR^2}{4k}\sum_{k=1}^{\infty}\frac8{\lambda_k^3J_1(\lambda_k)}J_0\left(\frac{\lambda_kr}R\right)\left(1-e^{\frac{\lambda_kk}{\rho C_p}t}\right)$$
