I am attempting to solve the 1D Heat Equation in cylindrical coordinates using separation of variables. Overall, I feel like I understand everything perfectly up until the part where I have to apply boundary conditions to find the values of my integration constants. I was hoping someone could clarify how to do this in more detail, and where to go from here.
The 1D Heat Equation is defined as: $$\frac{dT}{dt} = \frac{k}{\rho C_p}(\frac{1}{r}\frac{dT}{dr}+\frac{d^2T}{dr^2}).$$ Let $$A =\frac{k}{\rho C_p}.$$
Using separation of variables, I managed to work out the following equations:
$$T(r,t) = X(r)\,Y(t)$$ $$\frac{Y'(t)}{Y(t)} + A\,m^2 = 0 $$ $$rX'(r)+r^2X''(r) + r^2m^2X(r) =0$$
These equations appear to give fairly straightforward solutions of: $$Y(t) = C_1e^{-Am^2t} $$ $$X(r) = C_2J_0(mr)+C_3Y_0(mr)$$
I understand everything up to here. However, I am lost as to how to properly apply the following boundary conditions to find the values of both the constants and $m$:
$$T(r,0) = Ti$$ $$T(R,t) = T_\text{surface}$$ $$\frac{dT(0,t)}{dr} = 0$$
My understanding is that somehow the equations for $X$ and $Y$ turn into infinite series, and you have to solve for an infinite number of constants and $m$ values (which I assume you can limit to some reasonable value like 100), but I have no idea how to actually do that - efficiently or otherwise. I also don't understand how to get from my equations for $X$ and $Y$ to the infinite series that I've seen elsewhere: it seems arbitrary to just "add" a summation term.
Most textbooks I've seen describe the portion above in great detail, but gloss over everything after it and just kind of say "Solve it now - good luck!". I was hoping to use this analytical approach as a faster/more accurate alternative to the finite difference method, but using root-finding methods on an infinite series sounds neither faster or more accurate than finite difference intuitively, since you will have so many constants to solve for that each will have their own residuals.
Firstly, to do separation of variables, you need to have homogeneous boundary conditions. You can get this by defining a new variable $$ u(r,t) = T(r,t) - T_\text{surface} $$ So that $u(r,t)$ will satisfy the same PDE, but with condition $u(R,t) = 0$, and initial condition $u(r,0) = T_i - T_\text{surface}$.
Now you apply separation of variables by assuming $u_n(r,t) = X_n(r)Y_n(t)$ (I subscript with an $n$ here to indicate we're going to find an infinite number of separable solutions, then argue any linear combination also satisfies the given PDE and boundary conditions, which is the 'adding a summation' bit). You can satisfy the boundary condition $u(R,t) = 0$ by requiring $X_n(R) = 0$ (and the inner boundary condition, which is really required to make sure the temperature remains bounded, with $X_n'(0) = 0$.
Solving the problem for $X$ requires a bit of knowledge about Bessel functions. The condition at $r=0$ says that $c_2=0$, as $Y_0$ has a logarithmic singularity at zero. The condition at $r=R$ specifies that $m = m_n$ is the $n$th root of $$ J_0(mR) = 0. $$ This is the bit you would have to work out numerically.
After this point, and having found $Y_n$, you argue that any linear combination of $u_n(r,t)$ must satisfy the PDE and boundary conditions (this requires that the PDE and boundary conditions are homogeneous). So you have $$ u(r,t) = \sum_{n=1}^\infty A_n\mathrm e^{-m_n^2At} J_0(m_nr) $$ The last step is applying the initial condition. To do this you need to know how to expand the initial condition in terms of the Bessel functions. This is a generalisation of the idea of Fourier series. You should look up Sturm-Liouville theory to understand how this works. The Bessel functions are orthogonal under the weight function $r$, so $$ A_n = \frac{\langle f(r), J_0(m_nr)\rangle}{\langle J_0(m_nr), J_0(m_nr)\rangle} = \frac{\int_0^R rf(r)J_0(m_nr) \,\mathrm dr}{\int_0^R rJ_0(m_nr)^2 \,\mathrm dr}, $$ where in this case $f(r) = T_i - T_\text{surface}$. This part you may also have to do numerically. Of course you can add $T_\text{surface}$ on at the end to get back to $T(r,t)$.
Regarding the value of an analytic vs numerical solution: the best thing about an analytical solution is that you can analyse it. As an example, even if you don't know the roots $m_n$ exactly, the analytic solution tells you that for large time the solution will approach the term with the largest (least negative value) of $m$, as the higher terms decay more quickly. Finite differences wouldn't tell you that.