You have to view Green's functions as inverse operators. In your case, you want $\Phi$, so the operator you wish to invert is the RHS:
$$
\Delta-\epsilon M
$$
so the Green's function you are looking for is:
$$
(\Delta_x-\epsilon M(x))G(x,y)=\delta(x-y)
$$
Note that since $M$ breaks the translation invariance, $G$ takes two spatial arguments, and does not depend only on the difference. $\Phi$ can therefore be expressed as:
$$
\Phi(x) = \int d^Dy G(x,y)j(y)
$$
so that:
$$
(\Delta-\epsilon M)\Phi = j
$$
in particular, you can take $j=M$.
Since you are assuming that $\epsilon$ is a small parameter, you can write $G$ as a perturbative series in $\epsilon$. This is the Born approximation if you truncate at leading order and interpret $M$ as a potential energy and $\Phi$ a first quantisation wave function.
Formally, using the unperturbed case:
$$
\begin{align}
G_0 &= \Delta^{-1} \\
&= -\frac{1}{4\pi |x-y|}
\end{align}
$$
you would write using linear algebra:
$$
\begin{align}
G &= \frac{1}{\Delta-\epsilon M} \\
&= \frac{1}{\Delta}\frac{1}{1-\epsilon \Delta^{-1}M} \\
&= \frac{1}{1-\epsilon \Delta^{-1}M}\frac{1}{\Delta} \\
&= (1+\epsilon \Delta^{-1}M)\frac{1}{\Delta} \\
&= G_0 +\epsilon G_0 MG_0 \\
G(x,y) &= G_0(x,y)+\epsilon\int d^D z G_0(x,z)M(z)G_0(z,y) \\
&= -\frac{1}{4\pi |x-y|}+\frac{\epsilon}{(4\pi)^2}\int d^D z \frac{M(z)}{|x-z||y-z|}
\end{align}
$$
with the associate diagrammatic representation.
Hope this helps.
