Laplace PDE with impulse function I want to solve the PDE: $$u_{xx}+u_{yy}=-\delta(x-x_0)\delta(y-y_0)$$ 
in the rectangle: $\quad R=\{(x,y):0\leq x\leq a\quad0\leq y\leq b\},\quad (x_0,y_0)\in R$
with the boundary conditions: $$u_x(0,y)=u_x(a,y)=0$$
$$u(x,0)=u(x,b)=0$$
Which method should I use when there is an impulse function (in this case Dirac's delta)?
 A: When Dirac delta functions appear in a PDE, this PDE shall be understood in the weak sense. That is, one needs to define a solution space $\Omega_1$ to which the solution of this PDE belongs, and some test space $\Omega_2$ with which the Dirac delta functions could be interpreted appropriately. In particular, the choice with $\Omega_1=\Omega_2$ is among the most popular ones.
Provided the boundary conditions
\begin{align}
u_x|_{x=0}=u_x|_{x=a}&=0,\\
u|_{y=0}=u|_{y=b}&=0,
\end{align}
it is appropriate to define
$$
\Omega_1=\text{span}\left\{\cos\frac{m\pi x}{a}\sin\frac{n\pi y}{b}:m,n\in\mathbb{N}\right\},
$$
because any continuous function in this space satisfies the assigned boundary conditions. Further, as we expect $u\in\Omega_1$, we have
$$
u(x,y)=\sum_{m,n=0}^{\infty}\hat{u}_{mn}\cos\frac{m\pi x}{a}\sin\frac{n\pi y}{b}.
$$
The way that the Dirac delta function is then interpreted in this way: $\forall\,v\in\Omega_2$, $v$ yields
$$
\left(-\Delta u,v\right)=\left(\delta(\cdot-x_0)\delta(\cdot-y_0),v\right),
$$
where $\left(\cdot,\cdot\right)$ is a bilinear functional on $\Omega_1\times\Omega_2$, i.e.,
$$
\left(f,g\right)=\int_0^a{\rm d}x\int_0^bf(x,y)g(x,y){\rm d}y,
$$
and $v$ is always chosen to be smooth or even analytic. Obviously, this bilinear functional is no more than an integration.
Now, take $\Omega_2=\Omega_1$, and set
$$
v=\cos\frac{p\pi x}{a}\sin\frac{q\pi y}{b}
$$
with $p,q\in\mathbb{N}$ being arbitrarily chosen. With this choice, the weak formula
$$
\left(-\Delta u,v\right)=\left(\delta(\cdot-x_0)\delta(\cdot-y_0),v\right)
$$
implies
$$
\frac{\pi^2ab}{4}\left(\frac{p^2}{a^2}+\frac{q^2}{b^2}\right)\hat{u}_{pq}=\cos\frac{p\pi x_0}{a}\sin\frac{q\pi y_0}{b}.
$$
Therefore, we have
$$
\hat{u}_{mn}=\frac{4}{\pi^2ab}\left(\frac{m^2}{a^2}+\frac{n^2}{b^2}\right)^{-1}\cos\frac{m\pi x_0}{a}\sin\frac{n\pi y_0}{b}.
$$
This gives
$$
u(x,y)=\sum_{m,n=0}^{\infty}\left[\frac{4}{\pi^2ab}\left(\frac{m^2}{a^2}+\frac{n^2}{b^2}\right)^{-1}\cos\frac{m\pi x_0}{a}\sin\frac{n\pi y_0}{b}\right]\cos\frac{m\pi x}{a}\sin\frac{n\pi y}{b}.
$$
I would like to thank @Dylan for pointing out my arithmetic mistakes here.
