Modified Bessel equation with Gaussian source I'm studying a system in which a neutral particle beam with a Gaussian profile of known width $w$ enters a plasma. The beam produces secondary particles, which diffuse and are ionised through a first-order reaction (volume rate proportional to density $n$ of secondaries). The equation governing the density profile of secondary particles away from a line source is the modified Bessel equation,
$$r\frac{d}{dr}\left(r\frac{dn}{dr}\right) - \lambda^2r^2n = 0 \; , $$
which is solved by modified Bessel functions of the second kind, $K_0(\lambda r)$. My problem comes when I try to consider the profile inside the source beam; the RHS of the above equation becomes a Gaussian, $A\exp(-r^2/w^2)$, with $A$ a known constant.
At first, I thought I'd be able to convolve the solution to the Bessel equation with the Gaussian to get the overall profile, but when it turned out that the Bessel function has a singularity at the origin I realised this would not be possible. I haven't been able to find any suitable solutions online, so I'd be very grateful for any insight the SE community can provide.
 A: The question is about a mathematical method to solve a second order ODE describing a physical situation. The mathematical problem is of a standard nature: given a homogeneous, linear, second-order ODE, what's the solution of the inhomogeneous equation? 
The general solution of the homogeneous modified Bessel equation (with index 0) is (using the notation of the question)
$$n(r) = c_1 I_0 (\lambda r) + c_2 K_0 (\lambda r) \; .$$
Note that it has two independent ("fundamental") solutions as is appropriate for a second order ODE. $c_1$ and $c_2$ are fixed by boundary conditions.
The general solution of the inhomogeneous equation is found in most texts on ODEs, or online here. Denoting the inhomogeneity by $g(r)$ (which in the question is a Gaussian), the solution becomes 
$$n(r) = - I_0 (\lambda r) \int du \, \frac{K_0(u)g(u)}{W(I_0, K_0)} + K_0 (\lambda r) \int du \, \frac{I_0(u)g(u)}{W(I_0, K_0)} \; , $$
where $u = \lambda r$ and $W$ denotes the Wronskian,
$$W(I_0, K_0) = I_0 K_0' - K_0 I_0' \ne 0 \; ,$$
the prime denoting a derivative with respect to the argument $(\lambda r)$. The integration boundaries are to be chosen such that $n(r)$ obeys the appropriate boundary conditions.
I'd be surprised if the integrals could be done analytically, but numerically they should be feasible.
