Solution 1:

If you do not assume radial symmetry, and you let $T_{eq}=R(r)\Theta(\theta)$, then by using the separation of variables, you will arrive at a Cauchy-Euler equation for $R$.

I think because $T_{eq}$ satisfies Laplace's equation rather than the heat equation, it is wrong to assume radial symmetry.

So use separation of variables, to obtain:

$\frac{1}{r}R'\Theta+R''\Theta+\frac{1}{r^{2}}R\Theta'',\,\,\,$ now divide by $R\Theta/r^2$

$\frac{rR'+r^2R''}{R}+\frac{\Theta''}{\Theta}=0$

Solution 2:

The problem that the OP is trying to solve can be made non-dimensional by defining the following quantities:

$$\phi = \frac{T - T_0}{T_1 - T_0}, \quad \xi = \frac{r}{a}, \quad \tau = \frac{\kappa t}{a^2}, $$ hence obtaining (prove it):

$$ \frac{\partial \phi}{\partial \tau} = \frac{1}{\xi} \frac{\partial }{ \partial \xi} \left( \xi \frac{\partial \phi}{\partial \xi} \right), \quad 0 < \xi < 1, \quad \tau > 0,$$ where we have assumed that the temperature distribution is axisymmetric, since $\partial_\theta = 0$ everywhere (see comments). This equation is to be solved with the following boundary and initial conditions:

$$\begin{align} |\phi(0,\tau ) | & \leq \infty, \\ \phi(1,\tau) & = 1, \\ \phi(\xi,0) & = 0. \end{align} $$ Since the boundary conditions are not homogenous, we can make use of the superposition principle and set $\phi = u + v$, where $u$ satisfies homogenous boundary conditions and $v$ satisfies $|v(0,\tau)| < \infty$ and $v(1,\tau) = 1$. We can therefore set $v = 1$ as a possible solution for $v$. The problem for $u$ is then given by:

$$ \frac{\partial u}{\partial \tau} = \frac{1}{\xi} \frac{\partial }{ \partial \xi} \left( \xi \frac{\partial u}{\partial \xi} \right), \quad 0 < \xi < 1, \quad \tau > 0,$$

$$\begin{align} |u(0,\tau ) | & \leq \infty, \\ u(1,\tau) & = 0, \\ u(\xi,0) & = T(\xi,0) - v(\xi,0) = -1, \end{align} $$ where we have made use of $\partial_t v= \partial_\xi v= 0$. Since the PDE for $u$ is linear, we can use the method of separation of variables to write: $u(\xi, \tau) = P(\xi)Q(\tau)$, for $P, Q$ non-zero functions of their respective arguments. Introducing this decomposition into the PDE yields:

$$ P Q' = \frac{Q}{\xi} \frac{d}{d \xi} \left( \xi \frac{d P}{d\xi} \right),$$ which can be reduced to:

$$ \frac{Q'}{Q} = \frac{1}{\xi \, P} \frac{d}{d \xi} \left( \xi \frac{d P}{d\xi} \right) = \lambda \in \mathbb{R}.$$ This equations represents two problems: one for $P$, which is the interesting one, and one for $Q$ which we will not worry about. The problem for $P$ reads:

$$\frac{d}{d \xi} \left( \xi \frac{d P}{d\xi} \right) - \lambda \xi P = 0, \quad 0 < \xi < 1, $$ and $|P(0)|< \infty, \ P(1) = 0$ (see the definition of $u$ to derive this boundary conditions). This equations has non-trivial BCs-satisfying solution for $\lambda < 0$ if I remember well, so that $\lambda = - |\lambda|$, and then:

$$P(\xi) = C_1 \, J_0 \left(\sqrt{|\lambda |} \xi \right) + C_2 \, Y_0 \left(\sqrt{| \lambda | } \xi \right),$$ where $J_0$ and $Y_0$ are the Bessel functions of first and second kind of order zero, respectively. Since $Y_0$ is singular at $\xi = 0$ and we must have bounded solutions for both $P$ and $u$, the constant of integration $C_2$ must vanish. The condition $P(1) = 0$ leads to:

$$ 0 = C_1 \, J_0 \left(\sqrt{|\lambda |} \right), $$ which holds whether $C_1 = 0$ (which yields the trivial, and not desired, solution $P = 0$) or $\sqrt{|\lambda|} = \chi_n$ is a zero of $J_0$. Note that $J_0$ has infinitely many zeros for $\chi_n = 0$. This values of $\chi_n$ are called the eigenvalues of the problem and allows us to write the solution as:

$$P(\xi) = C_n J_0 (\chi_n \xi) = J_0 (\chi_n \xi), \quad n = \{1,2, \ldots \} $$ where I have made $C_n = 1$ for reasons you are about to see.

The Sturm-Liouville theory tells us that we can expand the solution as a linear combination of the eigenfunctions as follows:

$$u(\xi, \tau ) = \sum_{n=1}^{\infty} Q_n(\tau) P_n(\xi),$$ where $P_n(\xi ) = J_0 (\chi_n \xi)$ and $Q_n(\tau)$ are the Fourier coefficients of the expansion, which are to be determined introducing this into the original PDE for $u$:

$$ \sum_{n=1}^{\infty} Q'_n(\tau) P_n(\xi) = \sum_{n=1}^{\infty} Q_n(\tau) \underbrace{\frac{1}{\xi} \frac{d}{d \xi} \left( \frac{1}{\xi} \frac{d P_n}{d \xi} \right) }_{- |\lambda_n| }$$

Work in progress