I am trying to solve the heat equation in a disk of radius $R$ with Neumann (no flux) boundary conditions and Dirac distributed initial condition and radial symmetry. Starting from the backward Fokker-Planck equation for the probability $P(r,t)$, the problem is: $$\begin{matrix}\frac{\partial P}{\partial t}=D\left(\frac{\partial^2 P}{\partial r^2}+\frac{1}{r}\frac{\partial P}{\partial r}\right)\\ \frac{\partial P}{\partial r}\Big|_{r=R}=0 \end{matrix} $$ subject to initial condition where P is concentrated at the center $r=0$ which mathematically is $$ P(r,0)=\delta(r) $$
Based on separation of variables I managed to get that the general solution reads $$P(r,t)=\sum_{k=1}^\infty c_kJ_0(\lambda_kr)e^{-\lambda_k^2Dt}$$ Then I elucidated the eigenvalues $\lambda_k$ by plugging this expression into Neumann's boundary condition and obtained that $\lambda_k=\alpha_{1,k}/R$ where $\alpha_{1,k}$ is the k-th zero of $J_1(x)$. $$P(r,t)=\sum_{k=1}^\infty c_kJ_0(\alpha_{1,k}r/R)e^{-\alpha_{1,k}^2Dt/R^2}$$
Intermediate question: I expect that at long time, my solution converges to a uniform distribution over the disk, how can that be when each term contains $e^{-\alpha_{1,k}^2Dt/R^2}$ that will eventually vanish? Should there not be an additional constant in the solution? Where would it come from?
Ignoring this issue, I then tried to idenfity the coefficients $c_k$, I multiplied both sides of initial condition $rJ_0(\alpha_{1,m}r/R)$ and integrate over $R$ $$ \sum_{k=1}^\infty c_k \int_0^R rJ_0(\alpha_{1,k}r/R)J_0(\alpha_{1,m}r/R)dr=\int_0^R\delta(r)rJ_0(\alpha_{1,m}r/R)dr$$
Considering the left hand side: $\alpha_{1,k}$ not being a zero of $J_0(x)$, after some research I found that the left-hand side is a 'Lommel integral': For $\alpha_{1,m}\neq \alpha_{1,k}$ the integral is equal to zero (so we get rid of the infinite sum), but for $\alpha_{1,m}= \alpha_{1,k}$ I get $$ \int_0^R r[J_0(\alpha_{1,k}r/R)]^2dr= R^2/2[J_0(\alpha_{1,k})]^2$$
Here comes the problem: on the right-hand side the integral containing $\delta(r)r$ is equal to 0 from which I get $c_k=0$ that makes no sense.
I have been working on it for days and I can't see my mistake. Could someone help me please? Thank you very much.
A point source at the origin is
$$ \delta(x)\delta(y)=\frac{1}{2 \pi r}\delta(r) \neq \delta(r) $$
The delta picks up an inverse Jacobian under co-ordinate transforms. I presume your separation constant was $\lambda^2$. Notice that for $\lambda^2=0$ the solution is $P(r,t)=P_0+b \ln r$. The expression in your question is for $\lambda \neq 0$. The general solution that is finite at the origin for all $t$ is
$$ P(r,t)=P_0+\sum\limits_\lambda c_\lambda e^{-D\lambda^2 t}J_0(\lambda r) $$
From which we see that $P_0$ is the steady state constant. I'll set $R=1$ and $D=1$ for clarity. The Neumann BC gives (as you say) $\lambda_n=\alpha_n$, roots of $J_1$.
$$ P(r,t)=P_0+\sum\limits_{n>0} c_ne^{-\alpha_n^2 t}J_0(\alpha_n r) $$
Evaluating at $t=0$ for the initial condition
$$ P(r,0)=P_0+\sum\limits_{n>0} c_n J_0(\alpha_n r) = \frac{\delta(r)}{2\pi r} $$
Multiply by $2 \pi r J_0(\alpha_{m}r)$ and integrate from $0$ to $R$
$$ 2 \pi P_0 \int\limits_0^1 dr \ r J_0(\alpha_{m}r)+2 \pi \sum\limits_{n>0} c_n \int\limits_0^1 dr \ r J_0(\alpha_n r) J_0(\alpha_{m}r)=1 $$
The first integral in non-vanishing only if $\alpha_m=0$, the second vanishes unless $n=m$. When $m=0 \neq n$ we find $P_0=\frac{1}{\pi}$, and when $m=n \neq 0$ the integral over Bessel functions may be performed, leaving
$$ c_n=\left[ \pi J_0(\alpha_n) J_0(\alpha_n) \right]^{-1} \qquad ;\qquad n=1,2,3 \dots $$