You can use the Fourier Transform theory to seek for a solution. I think you are wrong in making your inverse transforms.
The first task is to transpose your problem into Fourier space. We shall use vector notation letting first $\vec x = (x,y)$, $\vec k=(k_x,k_y)$ the wave vector conjugated in 2D Fourier plane with spatial position $\vec x$ in direct plane. Let's write $\delta(\vec x-\vec x_0)$ the 2D delta distribution in direct plane. We define the direct and inverse Fourier transform of a distribution of
density $f(x)$ by,
$$\hat f(\vec k)=\int f(\vec x) \exp(-i \vec k.\vec x) d\vec x$$
and its inverse (normalization of the FT is inessential in what follows),
$$f(x)=\int \hat{f}(\vec k) \exp(i \vec k.\vec x) d\vec k$$
You seek for solving the Poisson equation,
$$\nabla^2 v(\vec x) = -\frac{q}{\epsilon} \delta(\vec x-\vec x_0)$$
It is very easy to show that the action of the Laplacian operator $\nabla^2 v(\vec x)$ translates on Fourier plane into a multiplication of $\hat v(\vec k)$ with the opposite of the squared modulus of the wave vector $-k^2$. Transforming the RHS of the equation is also straightforward when using the basic properties of the Dirac function. Hence your equation translates on Fourier plane into,
$$\hat v(\vec k)=\frac{q}{\epsilon} \frac{1}{k^2} \exp(-i \vec k.\vec x_0)$$
Take the inverse Fourier transform of the preceding relation side by side,
$$ v(\vec x, \vec x_0) = \frac{q}{\epsilon} \int \frac{1}{k^2} \exp[-i \vec k.(\vec x-\vec x_0)] d\vec k$$
Make use of a system of polar coordinates on the Fourier plane writing the wave vector $\vec k = k \hat{e}_k(\theta)$ with $\hat e_k(\theta)$ the local radial unit vector of cartesian coordinates $(\cos\theta,\sin\theta)$ on the Fourier plane, where $\theta$ is the polar angle. Using this new system of coordinates, the inverse Fourier transforms into,
$$ v(\vec x, \vec x_0) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp(-i k \hat e_k(\theta).(\vec x-\vec x_0) d\theta$$
Similarly, make use of a system of polar coordinates on the spatial plane to write down the displacement $(\vec x-\vec x_0)$ as $(\vec x-\vec x_0)=r \hat e_r(\phi)$ with unit radial vector $\hat e_r(\phi)=(\cos \phi,\sin\phi)$ where $\phi$ is the polar angle on the spatial plane. This yields,
$$v(\vec x,\vec x_0) \equiv v(r=|\vec x-\vec x_0|,\phi) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp[-i kr (\hat e_k(\theta).\hat e_r(\phi))] \ d\theta$$
The scalar product $\hat e_k(\theta).\hat e_r(\phi)$ is simply $\cos(\theta-\phi)$. It is always possible to choose the referent cartesian frame of the polar system such that $\phi=0$. Hence the transform is isotropic (does not depend upon $\phi$) reducing to,
$$ v(r) = \frac{q}{\epsilon} \int_{0}^{\infty} dk \int_{0}^{2\pi} \frac{1}{k} \exp[-i kr \cos\theta] \ d\theta$$
Integration over polar angle $\theta$ is straightforward it gives the Bessel function of the first kind $2\pi J_0(kr)$. As a conclusion you obtain an integral representation of your solution.
$$v(r)=\frac{2\pi q}{\epsilon} \int_{0}^{\infty} \frac{ J_0(kr)}{kr} d(kr)$$
Note that this integral representation is divergent. Introducing a strictly positive lower bound $a>0$ in the low part of the wavenumber spectrum allows to recover a definite solution. Further integrating over new variable $s=kr$,
$$v(r,a)=\frac{2\pi q}{\epsilon} \int_{ar}^{\infty} \frac{ J_0(s)}{s} ds$$
For instance mathematica/wolfram gives a function which is equivalent at the vicinity of $a=0$ to,
$$v(r,a) \sim -\frac{2\pi q}{\epsilon} [\log(\frac{1}{2} ar) +\gamma]$$
where $\gamma$ is the Euler-Mascheroni constant.
