I don't think that it is solvable for an arbitrary function $\epsilon(z)$. Here is how I would try to solve it in general case.
First of all, let us write it explicitly as
$$
-\epsilon(z)\left[\partial_x^2\phi(\mathbf{r},z) + \partial_y^2\phi(\mathbf{r},z)\right] -\partial_z\left[\epsilon(z)\partial_z\phi(\mathbf{r},z)\right] =
4\pi\delta(\mathbf{r} - \mathbf{r}_0)\delta(z-z_0).$$
Since the equation is homogeneous in transversal direction, we can use Fourier transform:
\begin{array}
a\phi(\mathbf{r},z) = \int\frac{dk_xdk_y}{(2\pi)^2}e^{-i\mathbf{k}\mathbf{r}}\tilde{\phi}(\mathbf{k},z),\\
\delta(\mathbf{r} - \mathbf{r}_0) =
\int\frac{dk_xdk_y}{(2\pi)^2}e^{-i\mathbf{k}(\mathbf{r}-\mathbf{r}_0)}.
\end{array}
Plugging these into the original equation we obtain
$$
-\partial_z\left[\epsilon(z)\partial_z\tilde{\phi}(\mathbf{k},z)\right]
+ \mathbf{k}^2\epsilon(z)\tilde{\phi}(\mathbf{k},z) =
4\pi\delta(z-z_0)e^{i\mathbf{k}\mathbf{r}_0}.
$$
This is a second order equation of type
$$
\frac{d}{dx}\left[p(x)\frac{d}{dx}y(x)\right] - k^2p(x)y(x) = 0.
$$
In general case this equation is not solvable, but it has known solutions for many types of function $p(x)$, since it is a Sturm-Liouville equation with zero eignevalue.
Assuming that we know two linearly independent solutions of this equation, $f_k(x)$ and $g_k(x)$, such that $f_k(x)\rightarrow 0$ as $x\rightarrow -\infty$ and $g_k(x)\rightarrow 0$ as $x\rightarrow +\infty$, we can write the solution of our equation of interest as
$$
\tilde{\phi}(\mathbf{k},z) =
\begin{cases}
Af_k(z), \,\,\,\, z<z_0,\\
Bg_k(z), \,\,\,\, z>z_0.\
\end{cases}$$
The constants $A$ and $B$ can be obtained from the boundary conditions (the second of which is obtained by integrating the equation over an infinitesimal interval $[z_0-\eta, z_0 +\eta]$:
\begin{array}
\tilde{\phi}(\mathbf{k},z_0 - \eta) = \tilde{\phi}(\mathbf{k},z_0 + \eta),\\
\epsilon(z_0)\partial_z\tilde{\phi}(\mathbf{k},z_0-\eta) - \epsilon(z_0)\partial_z\tilde{\phi}(\mathbf{k},z_0+\eta) = 4\pi e^{i\mathbf{k}\mathbf{r}_0},
\end{array}
that is
\begin{array}
Af_k(z_0) = Bg_k(z_0),\\
A\epsilon(z_0)\partial_z f_k(z_0) -B\epsilon(z_0)\partial_z g_k(z_0) = 4\pi e^{i\mathbf{k}\mathbf{r}_0}.
\end{array}
With known $A$ and $B$ we are know in the position to reassemble the solution and calculate the Fourier transform to get $\phi(\mathbf{r},z)$. It is certainly a hard way to take for the case when $\epsilon(z)$ is a constant, but for specific shapes of $\epsilon(z)$ it may yield a solution in terms of relevant special functions.
Update: solution for a charged plane
An interesting solvable case is a plane of charge located at $z=z_0$, in which case the principal equation takes form:
$$-\nabla\left[\epsilon(z)\nabla\phi(\mathbf{r},z)\right] = 4\pi\sigma\delta(z-z_0),$$
where $\sigma$ is the surface charge density.
Due to symmetry in $xy$-plane the solution depends only on $z$, i.e. $\phi(\mathbf{r},z) = \phi(z)$, and the equation can be written as
$$-\frac{d}{dz}\left[\epsilon(z)\frac{d}{dz}\phi(z)\right] = 4\pi\sigma\delta(z-z_0).$$
The boundary conditions at $z=z_0$ include continuity of the potential, $\phi(z_0-\eta) = \phi(z_0+\eta)$, and the boundary condition for the electric field that can be obtained by integrating the equation withing infinitesimally small region around $z_0$:
$$-\int_{z_0-\eta}^{z_0+\eta}dz\frac{d}{dz}\left[\epsilon(z)\frac{d}{dz}\phi(z)\right] = \epsilon(z_0)\left[\phi '(z_0-\eta) - \phi '(z_0+\eta)\right] = \int_{z_0-\eta}^{z_0+\eta}dz\frac{d}{dz} 4\pi\sigma\delta(z-z_0) = 4\pi\sigma.$$
The solutions on both sides of the charged plane are:
\begin{array}
\phi\phi(z) = \begin{cases}
B + A\int_z^{z_0}\frac{dz'}{\epsilon(z')}, z < z_0,\\
D + C\int_{z_0}^z\frac{dz'}{\epsilon(z')}, z > z_0.
\end{cases}
\end{array}
Imposing the boundary conditions we obtain:
\begin{array}
BB = D = \phi_0,\\
-A = C + 4\pi\sigma = -E_+\epsilon(z_0) + 4\pi\sigma,
\end{array}
where we defined the potential at point $z=z_0$ and the electric field immediately to the right from the charged plane, $E_+ = -C/\epsilon(z_0)$.
The solution is thus
$$
\phi(z) = \begin{cases}
\phi_0 - \left[-E_+\epsilon(z_0)+ 4\pi\sigma\right]\int_z^{z_0}\frac{dz'}{\epsilon(z')}, z < z_0,\\
\phi_0 - E_+\epsilon(z_0)\int_{z_0}^z\frac{dz'}{\epsilon(z')}, z > z_0.
\end{cases}.
$$
The electric field is given by
$$
E(z) = -\frac{d}{dz}\phi(z) = \begin{cases}
\left[E_+\epsilon(z_0)- 4\pi\sigma\right]\frac{1}{\epsilon(z)}, z < z_0,\\
E_+\frac{\epsilon(z_0)}{\epsilon(z)}, z > z_0.
\end{cases}
$$
The known case of a charged plane is vacuum is obtained by setting $\epsilon(z)=1$, and assuming that there is no external electric field applied, so that we can assume by symmetry that the fields to the left and to the right of the charged plane have the same magnitude: $E_+=2\pi\sigma$.
Let us consider a special case with
$$\epsilon(z) = \begin{cases} \epsilon, z<0\\1, z>0\end{cases},$$
where we assume that $z<z_0$.
The electric field than can be written as
$$
E(z) = \begin{cases}
\frac{E_+- 4\pi\sigma}{\epsilon(z)}, z < z_0,\\
\frac{E_+}{\epsilon(z)}, z > z_0.
\end{cases} =
\begin{cases}
\frac{E_+- 4\pi\sigma}{\epsilon}, z < 0,\\
E_+- 4\pi\sigma, 0 < z < z_0,\\
E_+, z > z_0.
\end{cases}
$$
This field can be though of as created by two charge planes: the one at $z=z_0$ and the image plane at $z=z_0$ with the effective charge corresponding to the jump of the field at $z=0$:
$$4\pi\sigma_{eff} = E_+- 4\pi\sigma - \frac{E_+- 4\pi\sigma}{\epsilon} =
\frac{\epsilon-1}{\epsilon}\left[E_+- 4\pi\sigma\right].
$$
We can fix constant $E_+$ by demanding, as for a charged plane in vacuum, that the electric fields at $z=\pm\infty$ have the same magnitude, i.e. $(E_+-4\pi\sigma)/\epsilon=E_+$. We then obtain
$$E_+ = \frac{4\pi\sigma}{1+\epsilon}, \sigma_{eff} = \sigma\frac{1-\epsilon}{1+\epsilon}.$$
Remark
Note that the displacement vector $\mathbf{D}=\epsilon \mathbf{E}$ is determined only by the distribution of the free charges. It therefore would be tempting to take the known solution for a point charge to equation $\nabla\cdot \mathbf{D} = 4\pi\delta(\mathbf{r}-\mathbf{r}_0)\delta(z-z_0)$ and then obtain the electric field as $$\mathbf{E} = \frac{\mathbf{D}}{\epsilon(z)}.$$
This solution seems to be at odds with the insolubility of the potential equation stated above, as well as with the exactly solvable case of a sharp dielectric boundary $$\epsilon(z) =\begin{cases}\epsilon, z<0,\\ 1, z>0\end{cases}.$$
The resolution of this seeming paradox is in the fact that the (static) electric field should satisfy also the equation $$\nabla\times\mathbf{E}=0,$$
to which the suggested above simple solution does not satisfy! Note also that only vector field with zero curl can be represented as a gradient of a potential.
