Tough second order differential equation I can't figure out this diff equation (in cylindrical coordinate). How can I solve it ?
Any comments appreciated
$$ 
\frac{1}{r}\frac{d}{dr}(r\frac{dE}{dr})+\frac{d^2E}{dz^2}+(\epsilon_0 k_0^2-\frac{1}{r^2}k_x^2R^2)E=\epsilon k_0^2E_0
$$
$$
r_1<r<r_2 ;0<z<d   
$$
 $(r_1, r_2, d \in \Bbb R )$
$$E(r,0)=0; E(r,d)=0$$
$$E(r_1,z)=0; E(r_2,z)=0 $$
 A: Aside from the boundary conditions, the PDE and the inhomogeneous term doesn't depend on $z$.
To construct an ansatz which satisfy the boundary conditions at $z = 0$ and $d$, we 
rewrite $E$ as a sum of standing waves in the $z$-direction. More precisely,
$$E(r,z) = \sum_{n=1}^\infty E_n(r) \sin\left(\frac{n\pi z}{d}\right)\tag{*1}$$
For $z \in (0,d)$, the inhomogeneous term on the RHS can be rewritten as
$$\epsilon k_0^2 E_0 = \frac{4 \epsilon k_0^2 E_0}{\pi}\sum_{\ell=0}^\infty\frac{\sin\left((2\ell+1)\frac{\pi z}{d}\right)}{2\ell+1}\tag{*2}$$
This suggests in the ansatz $(*1)$, we only need to keep the terms for odd $n$.
From now on, we will assume $n = 2\ell + 1$ is an odd positive integer.
Throw $(*1)$ and $(*2)$ into the PDE, we obtain following ODE for $E_n(r)$ 
$$\left( r^2\frac{d^2}{dr^2} + r\frac{d}{dr} + \left( \epsilon_0 k_0^2 - \frac{n^2\pi^2}{d^2}\right) r^2 - k_x^2 R^2 \right) E_n(r) = \frac{4\epsilon k_0^2 E_0}{\pi n} r^2$$
When $\displaystyle\;n < \sqrt{\epsilon_0} \frac{k_0 d}{\pi}\;$, introduce variables $x, \nu$ and function $\phi_n(x)$:
$$x = r \sqrt{\epsilon_0 k_0^2 - \frac{n^2\pi^2}{d^2}},\quad
\nu = k_x R,\quad\text{ and }\quad
E_n(r) = \frac{4\epsilon k_0^2 E_0}{\pi n\left( \epsilon_0 k_0^2 - \frac{n^2\pi^2}{d^2}\right)} \phi_n(x)$$
we reduce the ODE for $E_n$ to an inhomogeneous 
Bessel equation in $\phi_n(x)$.
$$\left( x^2\frac{d^2}{dx^2} + x\frac{d}{dx} + (x^2 - \nu^2)\right) \phi_n(x) = x^2$$
Its general solution has the form
$$\phi_n(x) = A J_\nu(x) + B Y_\nu(x) + 
\frac{\pi}{2}\left(
Y_\nu(x)\int x J_\nu(x) dx - J_\nu(x)\int x Y_\nu(x) dx
\right)$$
where $A, B$ are constants to be determined. You can use your boundary condition at $r = r_1$ and $r_2$ to pin down $A$ and $B$.
When $\displaystyle\;n > \sqrt{\epsilon_0} \frac{k_0 d}{\pi}\;$,
the coefficient in the $r \to x$ change of variable  becomes imaginary,
one need to replace above expression by corresponding modified Bessel functions.
