We note here that all of the ensuing analysis uses notation that is interpreted in the sense of Distributions or Generalized Functions.
With that note, let's begin by breaking the problem down into components, each of which is hopefully elementary.
STEP 1:
First, we know that the Fourier Transform of the Dirac Delta $\delta$ is
$$\begin{align}
\mathscr{F}\{\delta\}(k)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\delta(x)e^{ikx}\,dx\\\\
&=\frac{1}{\sqrt{2\pi}}
\end{align}$$
ASIDE:
This implies that the Inverse Fourier Transform of the constant function $1$ is the Dirac Delta $\delta(x)$. We can rewrite this relationship as
$$\bbox[5px,border:2px solid #C0A000]{\int_{-\infty}^{\infty}e^{-ikx}\,dk=2\pi\,\delta(x)} \tag 1$$
where again the notation in $(1)$ is interpreted as a distribution.
STEP 2:
Second, recall that the Fourier Transform has the property such that the Fourier Transform of the $n$'th order derivative, $D^nf$, of a function $f$, is given by
$$\bbox[5px,border:2px solid #C0A000]{\mathscr{F}\{D^nf\}(k)=(ik)^n\mathscr{F}\{f\}(k)} \tag 2$$
STEP 3:
Now, define the ramp function $r(x)$ as
$$r(x)=
\begin{cases}
x&,x\ge 0\\\\
0&,x<0
\end{cases}$$
Note that the second derivative of the ramp function is $D^2r=\delta$. Using $(2)$ reveals that
$$\begin{align}
\mathscr{F}\{r\}(k)&=\frac{1}{(ik)^2}\mathscr{F}\{D^2r\}(k)\\\\
&=\frac{1}{\sqrt{2\pi}}\frac{1}{(ik)^2}
\end{align} \tag 3$$
STEP 4:
Taking the inverse Fourier Transform of $(3)$ and multiplying by $2\pi$ yields
$$\bbox[5px,border:2px solid #C0A000]{2\pi\,r(x)=\int_{-\infty}^{\infty}\frac{-1}{k^2}e^{-ikx}\,dk }\tag 4$$
STEP 5:
Finally, using $(4)$ we find that
$$\bbox[5px,border:2px solid #C0A000]{\frac{-1}{2\pi}\frac{q}{\epsilon_0}\int_{-\infty}^{\infty}\frac{1}{k^2}e^{ikx}\,dk=\frac{q}{\epsilon_0}r(x)}$$
We can use the preceding analysis to solve the more general one-dimensional Poisson Equation
$$f''(x)=\rho(x)$$
Taking Fourier Transforms yields
$$\hat f(k)=\frac{-1}{k^2}\hat \rho(k)$$
whereupon inversion reveals that
$$\begin{align}
f(x)&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\frac{-1}{k^2}\hat \rho(k)e^{-ikx}\,dk\\\\
&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\frac{-1}{k^2} \left(\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\rho(x')e^{ikx'}\,dx'\right)e^{-ikx}\,dk\\\\
&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\rho(x')\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\frac{-1}{k^2}e^{-ik(x-x')}\,dk\,dx'\\\\
&=\int_{-\infty}^{\infty}\rho(x')r(x-x')\,dx'\\\\
&=\int_{-\infty}^{x}(x-x')\,\rho(x')\,dx'\\\\
\end{align}$$
Thus, the general solution is
$$\bbox[5px,border:2px solid #C0A000]{f(x)=\int_{-\infty}^{x}(x-x')\,\rho(x')\,dx'}$$
Separation of variables gives
$$
\frac{X''}{X} = \lambda = -\frac{Y''}{Y}.
$$
Because you want solutions that remain bounded in $x$ as $|x|\rightarrow\infty$, that dictates $\lambda = -\mu^2$ where $\mu$ is real. Otherwise you get exponential solutions that explode at one or both of $\pm\infty$. Because $y(1)=0$ needs to hold, then the solutions are
$$
X(x)=e^{i\mu x},\;\;Y(y)=\sinh(\mu(y-1)).
$$
The trial solution is an integral "sum" of such solutions
$$
u(x,y)=\int_{-\infty}^{\infty}c(\mu)e^{i\mu x}\sinh(\mu(y-1))d\mu.
$$
The coefficient $c(\mu)$ is determined by
$$
u(x,0) = -\int_{-\infty}^{\infty}c(\mu)e^{i\mu x}\sinh(\mu)d\mu.
$$
The function $u(x,0)$ is $1$ for $-2 \le x \le 2$ and is $0$ otherwise. So you want to find coefficients $c(\mu)$ such that
$$
\chi_{[-2,2]}(x) = -\int_{-\infty}^{\infty}c(\mu)e^{i\mu x}\sinh(\mu)d\mu.
$$
Multiplying by $e^{-is x}$, integrating and using the Fourier orthogonality trick,
$$
\frac{1}{2\pi}\int_{-2}^{2}e^{-is x}dx = -c(s)\sinh(s) \\
\frac{1}{\pi}\frac{e^{-i2s}-e^{i2s}}{-2is}= -c(s)\sinh(s) \\
\frac{1}{\pi}\frac{\sin(2s)}{s}=-c(s)\sinh(s) \\
c(s) = -\frac{1}{\pi}\frac{\sin(2s)}{s\sinh(s)}.
$$
Therefore, the solution $u(x,y)$ is given by
$$
u(x,y)=-\frac{1}{\pi}\int_{-\infty}^{\infty}e^{i\mu x}\frac{\sin(2\mu)}{\mu}\frac{\sinh(\mu(y-1))}{\sinh(\mu)}d\mu.
$$
I may be off by a negative. I don't see how you're going to get a discrete sum out of that integral. You can easily check that the above is the correct solution, so far as it goes. Differentiating twice in $x$ gives the negative of differentiating twice in $y$. Clearly $u(x,1)=0$, and $u(x,0)$ isn't hard to check either. $\sinh(\mu)$ has zeros on the imaginary axis, and maybe they're using residues somehow, but I don't see how that would give arctan terms.
Best Answer
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.