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)?
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.