Difficulty in Solution of Poisson's equation using Fourier Transform I wanted to find the potential on an infinite plane due to a point charge located at some point $(x_0,y_0,z_0)$.
So I decided to solve the $2D$ Poisson's equation.
$$\nabla^2 V(x,y)=-\frac{q}{\epsilon}\delta (x-x_0)\delta (y-y_0)$$
Rather than the conventional Green's Function, my method of approach was Fourier transform.  Twice application of FT gave
$$V(k,l)=\frac{q}{2\pi \epsilon}\cdot \frac{\exp\left(ikx_0+ily_0\right)}{k^2+l^2}$$
Then I started doing the inverse transforms. After two inverse fourier transforms, I got the integral
$$V(x,y)=\frac{q}{4\pi \epsilon}\cdot \int_{-\infty}^{\infty} \frac{\exp\left(ik(x-x_0)+k\lvert y-y_0\rvert\right)}{k} dk$$
Now this integral on computation using complex analysis gives $\pi i$. And the potential becomes complex. Whereas it should have been logarithmic.
Can someone tell me where I am mistaken? And how I can better compute the integral?
Or maybe the problem lies in the first few steps??
 A: As Paul Garrett points out in the comments, solving a PDE with constant coefficients in the space of tempered distributions is the same as finding distributional inverses to polynomial functions. In this case, the solution involves finding a distribution $\gamma$ such that $\|x\|^2\cdot \gamma(x) = 1$. The natural guess $\gamma(x)=\frac{1}{\|x\|^2}$ does not work immediately because this is not locally integrable over the plane (where as in dimensions 3 and higher, this would work). We can only say that $\gamma(x)=\frac{1}{\|x\|^2}$ holds away from the origin.
Nevertheless, we can slightly modify the function to make work. We can add any linear combination of $\delta$ and $\frac{\partial}{\partial x_i}\delta$, since these are exactly the distributions with $\|x\|^2 u=0$.
Using a Taylor-expansion we can write: $\phi(x) = \phi(0) + \langle\nabla\phi(0),x\rangle + \tilde{\phi}(x)$ with $\tilde{\phi}\in o(\|x\|^2)$. And if we then remember than $\int_{\mathbb{R}^2} \frac{1}{\|x\|^2}\tilde{\phi}(x) dx$ is well-defined for Schwartz functions that vanish quadratically at the origin (because the integrand extends continuously to the origin), we can guess that a distributional inverse would be
$$\langle\gamma,\phi\rangle := \int_{\mathbb{R}^2} \frac{\phi-\phi(0)-\langle\nabla\phi(0),x\rangle}{\|x\|^2} dx$$
Now one can try to Fourier transform this distribution.
A: 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. 
A: Your equation for finding this potential is wrong since it is$$\nabla^2 V(x,y)=-\frac{q}{\epsilon}\delta(\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2})$$but i think solving this using the same relation of pointwise electrical charge potential is easier.$$V(x,y,z)=\dfrac{q}{4\pi\epsilon_0\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}}$$
