The Poisson's equation: $$ \nabla^2\phi=I\delta(\vec{r}-\vec{r_s}) $$ with the boundary condition $\phi(\infty)=0$, where $I$ is a constant and $\vec{r_s}$ is the source location, has the solution: $$ \phi(\vec{r})=\frac{I}{4\pi|\vec{r}-\vec{r_s}|} $$ I am trying to solve this equation in 2D by finite-difference method with grid sizes $dx=dy$. At each grid point not contaning the source, we have the equation: $$ \frac{1}{4}(\phi_{i,j-1}+\phi_{i,j+1}+\phi_{i-1,j}+\phi_{i+1,j})-\phi_{i,j}=0 \tag{*} $$ At the soure location, the RHS of (*) must include the source term.
Question: If the grid sizes are $dx=dy$, what is the value of the RHS of (*) at the source location?