Solving an IVBP for a metal rod dipped from cold to warm water I have to solve the diffusion equation for a solid rod, which is in a bath of 100 degrees. At $t=0$ it is moved to a second bath, where the water temperature is 0 degrees.
I prepare the diffusion equation:
\begin{equation}
u_t-\alpha\Delta=0.
\end{equation}
Boundary and initial conditions:
\begin{equation}
u(r,\theta,t)=0 \ \ \ \  0<\phi<\pi, t>0 \\
u(r,\phi,0)=100 \ \ \ \  0<r<R,  0<\phi<2\pi.
\end{equation}
Since the function $u$ is of three variables $u=R\Theta T$, I prepare the given partial differential equation with the differential operator  on a circle (the cross section of the rod):
\begin{equation}
T_tR\Theta-\alpha\bigg[R_{rr}\Theta T+\frac{1}{r}R_r\Theta T+\frac{1}{r^2}\Theta_{\theta\theta}RT\bigg]=0
\end{equation}
Divide by $R\Theta T$ and obtain
\begin{equation}
\frac{T_t}{T}-\alpha\bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=0
\end{equation}
Forming one  ODE and one PDE:
\begin{equation}
\frac{T_t}{T}=-\lambda^2 \\
\bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=\frac{\lambda_n^2}{\alpha}
\end{equation}
The first gives $T(t)=Ce^{-\lambda^2t}$, the second partial differential equation gives two ordinary differential equations:
\begin{equation}
r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}+\frac{\Theta_{\theta\theta}}{\Theta}=\frac{\lambda_n^2}{\alpha}r^2
\end{equation}
which can be further split into two ODEs which must equal some constant $\mu^2$, which must be negative to give Bessel solutions:
\begin{equation}
r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}-\frac{\lambda_n^2}{\alpha}r^2=\frac{\Theta_{\theta\theta}}{\Theta}=-\mu_k^2
\end{equation}
this gives the Bessel equation:
\begin{equation}
R_{rr}+\frac{1}{r}R_r+\bigg(\frac{\lambda_n^2}{\alpha}-\frac{\mu_k^2}{r^2}\bigg)R=0
\end{equation}
and the second-order ODE:
\begin{equation}
\Theta_{\theta\theta}+\mu_k^2\Theta=0
\end{equation}
The latter has the  solution $\Theta(\theta)=A\sin\mu_k\theta+B\cos\mu_k\theta$.
The former is the Bessel equation which gives the Bessel functions:
\begin{equation}
R(r)=\bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_n r}{\sqrt{\alpha}}\bigg)+bY_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)
\end{equation}
So if this procedure is right, the general solution would be (ignoring the Bessel function of the second kind:
\begin{equation}
u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\bigg(\big(A\sin\mu_k\theta+B\cos\mu_k\theta\big)\big(e^{-\lambda_n^2t})\bigg)
\end{equation}
Then we set the angular part equal to some constant, which is merged with the Bessel constant $a$ and get
\begin{equation}
u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\big(e^{-\lambda_n^2t})\bigg)
\end{equation}
Using ICs to find the constants, the zeros of the Bessel are given by $\alpha_{n,k}$ and since the core will still be at 100 degrees at t=0, while the surface will be 0, then $R(0)=100$ and $R(R)=0$, and with the radial function given by:
\begin{equation}
R(r)=aJ_{\lambda n}(i\sqrt{\mu}r)
\end{equation}
we solve:
\begin{equation}
\alpha=aJ_{\lambda n}(i\sqrt{\mu}R) 
\end{equation}
and
\begin{equation}
100=aJ_{\lambda n}(i\sqrt{\mu}\alpha) 
\end{equation}
which gives:
\begin{equation}
u(r,t)=\frac{\alpha J_{\lambda_n}(i\sqrt{\mu}r)}{J_{\lambda_n}\big(\frac{100R}{a\alpha}\big)}e^{-\lambda^2t}
\end{equation}
The plot of the time-dependent evolution looks like this (with x=r) and  x=0=R (surface) and t goes from 0 to 100. (ignore minus signs)

Does this entire procedure appear reasonable?
Thanks
 A: First let $u=X(\mathbf{x}) T(t)$. The equation reads $T_t X - \alpha \Delta(X) T = 0$. Take the cylindrical Laplacian and kill the $z$ and $\theta$ derivatives from symmetry considerations, so that $X(\mathbf{x})=R(\| \mathbf{x} \|)$; let $r=\| \mathbf{x} \|$.
Now you have
$$T_t R - \alpha T \left ( R_{rr} + r^{-1} R_r \right ) = 0.$$
Doing the usual separation of variables trick, you conclude that $\frac{T_t}{T}=\alpha \frac{R_{rr} + r^{-1} R_r}{R}=\lambda_n$ for the same constant $\lambda_n$. It follows that up to a constant factor that we'll sort out later, $T=e^{\lambda_n t}$.
Now in the radial part we have
$$R_{rr} + r^{-1} R_r - \frac{\lambda_n}{\alpha} R = 0.$$
This is actually an adjusted Bessel equation itself. Multiply through by $r^2$:
$$r^2 R_{rr} + r R_r + \left ( -\frac{\lambda_n}{\alpha} r^2-0^2 \right ) R = 0.$$
That $-\frac{\lambda_n}{\alpha}$ wouldn't be there in the ordinary Bessel equation. The fix comes because the first two terms are invariant under a linear change of variable, so we can replace $s=\sqrt{-\lambda_n/\alpha} r$ (dropping the $n$ as an abuse of notation) and then we have
$$s^2 R_{ss} + s R_s + \left ( s^2 - 0^2 \right ) R = 0.$$
This gives $R=c_{1,n} J_0(s) + c_{2,n} Y_0(s)$, Bessel functions of order zero. For nonsingularity at the origin we force $c_{2,n}=0$.
So now the question is about what $c_{1,n}$ is. The big catch is that in order for $J_0(\sqrt{-\lambda_n/\alpha} r)$ to have zeros, we must have $\sqrt{-\lambda_n/\alpha}$ real. (This also makes intuitive sense as it implies exponential decay.)
As a result you choose $\lambda_n$ so that $\sqrt{-\lambda_n/\alpha} R$ is a zero of $J_0$ (where $R$ is now the radius of the rod). We can pick the zero of interest to be $j_{0,n}$ which is where $n$ enters into the picture.
Now we solve that for the eigenvalue, and get $\lambda_n=-\frac{\alpha j_{0,n}^2}{R^2}$ which gives
$$u=\sum_n c_n J_0 \left ( j_{0,n} \frac{r}{R} \right ) e^{-\frac{\alpha j_{0,n}^2}{R^2} t}.$$
Now it remains to find the $c_n$ to match the initial profile.
