I am investigating diffusion of ink in water. I want to compare some experimental results to a solution $u=u(r,t)$ of polar diffusion equation $$D\cdot\left(\frac{\partial^2u}{\partial r^2}+\frac{1}{r} \frac{\partial u}{\partial r}\right)=\frac{\partial u}{\partial t},$$ (assuming isotropy, $\theta$ vanishes) with boundary conditions $$\begin{aligned}&u(r\le a, t=0)=u_0, \\ &u(r>a, t=0)=0.\end{aligned}$$ This post suggests a solution via separation of variables, however I cannot understand how it leads to the following solution, as shown in this article: $$u(r,t)=\frac{u_0}{2Dt}e^{\frac{-r^2}{4Dt}}\int_0^a{e^{\frac{-r^2}{4Dt}}I_0\left(\frac{rr'}{2Dt}\right)} r' dr',$$ Where $I_0(z)=J_0(iz)$ is the modified Bessel function of zeroth order.
Might this be a mistake? The repetition of $e^{\frac{-r^2}{4Dt}}$ inside and outside the integral looks suspicious, but this is the way it appears in the article. Can someone help to understand the derivation of the solution above, or give an easier to calculate solution? Much appreciated.
[Side note: If anyone is, by chance, familiar with a Matlab tool to solve this equation numerically (some parabolic equation solver perhaps - or, is solvepde helpful here?) it will be extremely helpful.]