Differential equation (inhomogeneous ) I have been trying to solve this equation for a while. Is there anyone who can help me to solve this ? Any comment appreciated.
$$\frac1r \frac{\partial}{\partial r}\left(r\frac{\partial E}{\partial r}\right) + \left(\frac{\partial^2E}{\partial z^2}\right) + \left(k_0^2-\frac{k^2}{r^2} \right) E(r,z) = f(x)$$
where $$r_1\leq r \leq r_2$$ $$0\leq z\leq d$$    $$E(r,0)=0$$ $$E(r,d)=0$$
 A: Assuming $x = z$, note that you can use separation of variables together with Sturm-Liouville theory as the boundary conditions in $z$ are homogeneous. First of all, let me consider the problem:
$$ \nabla^2 \phi + k_{\text{eff}}(r) \, \phi = 0, \quad (r,z) \in (r_1,r_2)\times(0,d), \quad \phi = \phi(r,z), $$ with homogenous BCs on $x$. Here, $\nabla^2$ stands for the laplacian operator in cylindrical axisymmetric coordinates, i.e., $\partial_\theta = 0$.
Make use of SV to set $\phi = P(r) Q(z),$ where $P$ and $Q$ are non-zero functions of their respective domain of definition. Introduce this in the PDE to come up with:
$$ - \mathcal{B}(P)/P = Q''/Q, \quad \mathcal{B} = \frac{1}{r} \frac{d}{dr}\left( r \frac{d}{dr} \right) + k_{\text{eff}}(r) .  $$
Since the equality holds for $P$ and $Q$ being functions of $r$ and $z$ respectively:
$$Q''/Q = \lambda \in \mathbb{R}, \quad Q(0) = Q(d) = 0,$$
where the BCs for $Q$ are derived from those of $\phi$ (remember they are homogenous in $z$). If we solve this 1st order ODE for the cases $\lambda = 0, \lambda < 0, \lambda > 0$ we finally find that the only non-zero solution occurs for $\lambda <0$, which is given by:
$$Q(z) = A \cos{\xi z} + B \sin{\xi z}, \quad \xi = \sqrt{|\lambda|}$$
Upon setting the boundary conditions we find:
$$ A = 0, \quad B \sin{\xi d} = 0, $$
the only solution for the last equation which permits $Q \neq 0$ is $\sin{\xi d} = 0$, that is to say:
$$ \xi d = n \pi \implies  \xi_n = \frac{n \pi}{d}, \quad n = 1, 2, \ldots$$
This values of $\xi_n$ are the eigenvalues of the problem in $z$. Therefore, by making use of SL theory we can expand the solution $E$ in terms of this eigenfunctions as follows:
$$ E(r,z) = \sum^\infty_{n=1} P_n(r) Q_n(z) = \sum^\infty_{n=1} P_n(r) \sin \xi_n z, $$
where $P_n$ are now the uknown coefficients of the Fourier expansion of $E$. Introduce now this result in your original PDE to have:
$$ \sum^\infty_{n=0} \left\{ \mathcal{B}(P_n) Q_n + P_n Q''_n \right\} = f(z) , $$
remember that $Q''_n = - \xi_2^n Q_n$ and $ \langle Q_i , Q_j\rangle_{z \in (0,d)} = \delta_{ij} \int_{(0,d)} Q_iQ_j \omega \, dz $, where $\omega = 1$ in your case and $\delta$ is the usual Kronecker delta. This allows us to obtain the final ODE for $P_n$ by multiplying the equation above by $Q_m \omega$ and integrating over the $z$ domain to have:
$$ \mathcal{B}(P_n) - \xi^2_n P_n = \displaystyle\frac{\displaystyle\int_0^d f(z) Q_n(z) \, dz  }{ \displaystyle\int_0^d Q^2_n(z) \, dz } = g_n(r).$$
The boundary conditions for $P$ are to be obtained from the definition of $E$ through the Fourier decomposition and using this very same trick.
The resulting equation for $P$ is a Bessel equation, namely:
$$ \frac{1}{r} \frac{d}{dr} \left( r \frac{dP_n}{dr} \right) + \left(  k_0^2 - \xi_n^2 - \frac{k^2}{r^2 } \right) P_n  = g_n,$$
whose homogenous solution is given by:
$$ P_n(r) = \alpha_n J_k(\mu_n r) + \beta_n Y_k(\mu_n r), $$ where $\alpha$ and $\beta$ are constants of integration, $\mu_n = k_0^2- \xi_n^2$ and $J_\nu$ and $Y_\nu$ are the Bessel functions of first and second kind of order $\nu$, respectively. 
I'm sure you can take it from here.
Cheers!
