I responded to a Quora question that asked about solving the electric field with radially decreasing permittivity. Specifically, solving a problem in the form, $\nabla^2 \textbf{E} + \omega^2 \mu_0 \varepsilon(0) e^{-kr} \textbf{E} =0$.

My answer was that the problem amounts to solving equations of the type,

$$(\nabla^2+f(r)){\boldsymbol{\mathrm{E}}}=0.$$

In Cartesian coordinates ($r^2=x^2+y^2+z^2$):

\begin{align}(\nabla^2+f(r))E_x&=0,\\

(\nabla^2+f(r))E_y&=0,\\

(\nabla^2+f(r))E_z&=0.\end{align}

If we want a propagating solution in the $z$-direction, we might as well set $E_x=E_y=0$. As for $E_z$, we get:

$$\frac{\partial^2}{\partial x^2}E_z+\frac{\partial^2}{\partial y^2}E_z+\frac{\partial^2}{\partial z^2}E_z+f(r)E_z=0.$$

We seek the solution in the form $E_z=\phi_x(x)\phi_y(y)\phi_z(z)$. Then, the equation becomes

$$\frac{\partial^2\phi_x}{\partial x^2}\phi_y\phi_z+\frac{\partial^2\phi_y}{\partial y^2}\phi_z\phi_x+\frac{\partial^2\phi_z}{\partial z^2}\phi_x\phi_y+f(r)\phi_x\phi_y\phi_z=0,$$

which translates into the following equations for $\phi_x$, $\phi_y$ and $\phi_z$:

\begin{align}\left[\frac{\partial^2}{\partial x^2}+\left(\frac{1}{\phi_y}\frac{\partial^2\phi_y}{\partial y^2}+\frac{1}{\phi_z}\frac{\partial^2\phi_z}{\partial z^2}+f(r)\right)\right]\phi_x&=0,\\

\left[\frac{\partial^2}{\partial y^2}+\left(\frac{1}{\phi_z}\frac{\partial^2\phi_z}{\partial z^2}+\frac{1}{\phi_x}\frac{\partial^2\phi_x}{\partial x^2}+f(r)\right)\right]\phi_y&=0,\\

\left[\frac{\partial^2}{\partial z^2}+\left(\frac{1}{\phi_x}\frac{\partial^2\phi_x}{\partial x^2}+\frac{1}{\phi_y}\frac{\partial^2\phi_y}{\partial y^2}+f(r)\right)\right]\phi_z&=0.\end{align}

In other words, the problem is reduced to solving an ODE in the form,

$$y''+f(x)y=0.$$

I do not know how to solve this type of equation in the general case. However, if $f(x)=A+e^{-kx}$, the solution is given in the form of Bessel functions:

$$f(x)=C_1J_\alpha\left(2\frac{e^{-kx/2}}{k}\right)+C_2Y_\alpha\left(2\frac{e^{-kx/2}}{k}\right),$$

where

$$\alpha=-2\frac{\sqrt{-A}}{k}.$$

This suggests to me that the solution to the problem in the general case would also be in the form of some generalized Bessel functions, i.e., a damped wave equation of some kind.