The question is from ISSP by Kittel and as follows:
(a)Find a solution of the London equation that has cylindrical symmetry and applies outside a line core. In cylindrical polar coordinates, we want a solution of $$ B-\lambda \nabla^2 B = 0 $$ that is singular at the origin and for which the total flux is the flux quantum: $$ 2 \pi \int_0^\infty d\rho\rho B(\rho) = \Phi_0$$ The equation is in fact valid only outside the normal core of radius $\xi$
(b) Show that the solution has the limits $$ B(\rho) \simeq (\Phi_0/2\pi\lambda^2)\ln(\lambda/\rho) , (\xi \ll \rho \ll \lambda) $$ $$ B(\rho) \simeq (\Phi_0/2\pi\lambda^2)(\pi\lambda/2\rho)^{1/2}\exp(-\rho/\lambda) , (\rho \gg \lambda) $$
I have proceeded to solve the equation with CAS software such as Maple and Mathematica; however, they give two arbitrary constants which are to be determined. Is the problem sufficiently defined to fix these arbitrary constants? The solutions given by these programs are: $$B(r)\to c_1 J_0\left(\frac{i r}{\lambda }\right)+c_2 Y_0\left(-\frac{i r}{\lambda }\right) $$ Furthermore, how can I describe the effect of core? I appreciate your help.
Edit: I found a solution manual but got confused by the method employed there in order to describe the effect of the core. Here is a screenshot of the page.

How can I implement 2-D Dirac delta function in cylindrical coordinates? By the way, I think one of the solutions of the homogenous equation is disregarded.
I won't do the whole question but I'll set it up for you and then hopefully everything will make sense.
Assume $B = R( \rho) \Theta ( \theta)$. Let's see how this changed the PDE, (you probably know how this goes so I'll skip to the punch line)
$$ B - \lambda \Delta B = 0 \implies - \frac{ \Theta''}{\Theta} = \rho^2 \frac{R''}{R} + \rho \frac{R'}{R} - \frac{\rho^2}{\lambda} $$
Notice that the right hand side only depends on $\rho$, and the left hand side depends on $\theta$, this can only be the case if this equals some constant $\mu \in \mathbb{R}$. This reduces the PDE to 2 ODE
$$ \Theta'' = - \mu \Theta \quad \& \quad \rho^2 R'' + \rho R' +( \mu + \frac{\rho^2}{\lambda} ) R = 0 $$
The $\Theta$ equation is our typical oscillation one and I just noticed we have cylindrical symmetry.... so $\Theta = const$..oh well, and the other is Bessel's equation. As you may have noticed this gives bessel type solutions! Now use the normilization you're given and see what you obtain. You'll be able to get all the constants, so part b shouldn't be too much of an issue now.
On second thought, you may want the following relation for the asymptotics :
$$J_n(x) = \frac{1}{\pi} \int_0^\pi \cos( n t - x \sin(t) ) dt $$