I've been working in 'Basic Partial Differential Equations' by Bleecker & Csordas' (Section 9.5, #5(a)) and I'm now having trouble with a this problem since the boundary conditions are left rather general. I'm not sure of how much more work I can do on the solution. The problem is as follows: Solve the heat equation in two dimensions \begin{align*} u_t &= k[u_{xx}+u_{yy}], \quad x^2+y^2< r_0, t\geq 0\\ u(r_0, \theta, t) &= 0, \quad -\pi \leq \theta \leq \pi\\ u(r,\theta,0) &= f(r,\theta) \end{align*}
After transforming our PDE into polar form we may use separation to find \begin{align*} \frac{T'}{kT}=\frac{R''\Theta+\frac{1}{r}R'\Theta+\frac{1}{r^2}R\Theta''}{RH}&=-\lambda^2\\ \implies \frac{r^2R''+rR'+\lambda^2r^2 R}{R} = -\frac{H''}{H} &= m^2 \end{align*}
where $-\lambda^2$ and $m^2$ are our constants of separation. Hence, \begin{align*} T'+k\lambda^2T&=0\\ r^2R''+rR'+(\lambda^2r^2-m^2)R&=0\\ H''+m^2H&=0 \end{align*}
We have general solutions to these ODEs \begin{align*} T(t) &= e^{-\lambda^2 k t} \\ R(r) &= c_1 J_m(\lambda r) + c_2 Y_m(\lambda r) \quad \text{(Bessel functions)} \\ \Theta(\theta) &= a_{m,q} \cos(m\theta) + b_{m,q}\sin(m\theta) \end{align*}
to form our general solution \begin{equation*} u(r,\theta,t) = \sum_q \sum_m u_{q,m}(r,\theta,t) = \sum_q \sum_m \left[a_{m,q} \cos(m\theta) + b_{m,q}\sin(m\theta)\right] \left[c_1 J_m(\lambda r) + c_2 Y_m(\lambda r) \right] e^{-\lambda^2 kt} \end{equation*} (i believe the indices range from $m = -\infty$ to $\infty$ and $q = 1$ to $\infty$, though I'm not quite sure why this is so)
Translating our boundary condition gives us \begin{equation*} u(r_0, \theta, t) = 0 \implies R(r_0)\Theta(\theta) = 0, \quad -\pi \leq \theta \leq pi \end{equation*}
Since it's not obvious to me that we know our problem is bounded (i.e. $|u(0)| < \infty$) to rule out second-kind Bessel functions $Y_m$, I don't know where to proceed in simplifying the solution any further, or what form our constants will take.
Thank you kindly!
Generally, the solution of a BVP must solve it at each point of domain. This in particular means that you should take the singular solution of your equation as a distribution. Neumann functions $Y_m$ don't solve the equation in polar coordinates in every point: they would only solve it if it had an additional (source) term proportional to Dirac delta with its singularity at the origin. See my answer at Physics.SE to a related question for extended discussion.