Finally figured out how to apply the method of images for this problem, it's similar as described previously in the same book (see page 83).
Consider the domain $D = \{(x_1,x_2): x_1 \in (-\infty,\infty), x_2 >0\}.$
We use the method of images to find the Green's function $g(x,\xi)$. Consider the point $(x_1,x_2)$ and it's image point $(x_1,-x_2)$. From (a), we add the function $$h = -\frac{1}{2\pi}\log \rho' = -\frac{1}{2\pi}\log \sqrt{(\xi - x_1)^2 + (\eta + x_2)^2}$$ to make $g = 0$ on the boundary. Since the image point is by definition not in the domain, $h$ is twice differentiable and satisfies Laplace's equation, $$\Delta h = \frac{\partial^2 h}{\partial \xi^2}+\frac{\partial^2 h}{\partial \eta^2}=0$$ for $(x_1,x_2)\in D$. Therefore, $g$ is defined by $$g = \frac{1}{2\pi}\log \rho + h = \frac{1}{2\pi}\log \rho - \frac{1}{2\pi}\log \rho' = \frac{1}{2\pi}\log{\frac{\rho}{\rho'}} = \frac{1}{4\pi}\log \frac{(\xi-x_1)^2 + (\eta - x_2)^2}{(\xi - x_1)^2 + (\eta+x_2)^2}.$$
On the boundary C, we have $\eta = 0$, therefore $g = 0$ and $$\nabla g = \frac{1}{\pi}\frac{x_2}{(x_1 - \xi)^2 + x_2^2}.$$ Inserting this solution of the two-dimensional problem into (1.1.5), we finally obtain $$u(x_1,x_2) = \frac1\pi \int_{-\infty}^{\infty}\frac{x_2 h(\xi)}{(x_1 - \xi)^2 + x_2^2}d\xi.$$
Note that (1.1.5) is given in the book as: $$u(x) = \int_\Omega g(x,\xi)f(\xi)d\xi - \int_\Gamma \frac{\partial g}{\partial n}h(\xi)dS_\xi.$$
We present an outline only for obtaining the solution to the posed problem for $u$ and provide an outline of a way for verifying the solution.
We proceed by using Green's Identity to write
$$\int_S (u(\vec r')\nabla'^2G(\vec r,\vec r')-G(\vec r,\vec r')\nabla'^2u(\vec r'))\,dS'=\oint_C (u(\vec r')\nabla' G(\vec x,\vec r')-G(\vec r,\vec r')\nabla' u(\vec r'))\cdot \hat n'\,d\ell' \tag 1$$
Using $\nabla' ^2G(\vec r,\vec r')=\delta(\vec r-\vec r')$ and $\nabla' ^2u(\vec r')=0$, the left-hand side of $(1)$ becomes
$$u(\vec r)=\int_S (u(\vec r')\nabla'^2G(\vec r,\vec r')-G(\vec r,\vec r')\nabla'^2u(\vec r'))\,dS' \tag 2$$
Using the boundary conditions,
$$\begin{align}
u(x,0)&=g(x)\\\\
G(x,0)&=0\\\\
\left.\frac{\partial u(x,y)}{\partial x}\right|_{x=0}&=h(y)\\\\
\left.\frac{\partial G(x,y;x',y')}{\partial x'}\right|_{x'=0}&=0\\\\
\end{align}$$
the integral on the right-hand side of $(1)$ can be written
$$\begin{align}
\oint_C (u(\vec r')\nabla' G(\vec x,\vec r')-G(\vec r,\vec r')\nabla' u(\vec r'))\cdot \hat n'\,d\ell'&=\int_0^\infty g(x')\left.\frac{\partial G(x,y;x',y')}{\partial y'}\right|_{y'=0}\,dx'\\\\
&-\int_0^\infty h(y')G(x,y;0,y')\,dy' \tag 3
\end{align}$$
Equating $(2)$ and $(3)$ yields the solution for $u(x,y)$ as
$$\bbox[5px,border:2px solid #C0A000]{u(x,y)=\int_0^\infty g(x')\left.\frac{\partial G(x,y;x',y')}{\partial y'}\right|_{y'=0}\,dx'-\int_0^\infty h(y')G(x,y;0,y')\,dy'} \tag 4$$
To check if $(4)$ is indeed a solution, we need to ensure that $u$ satisfies Laplace's equation for $x>0$, $y>0$ and the prescribed boundary conditions.
Since $G$ satisfies Laplace's equation for $(x,y)\ne (x',y')$, then $u$ as given in $(4)$ satisfies Laplace's equation.
Next, we examine $\lim_{y\to 0^+}u(x,y)$. Using heuristic formal analysis, we have for "small" $\delta >0$
$$\begin{align}
\lim_{y\to 0^+} u(x,y)&=\lim_{y\to 0^+}\int_0^\infty g(x')\left.\frac{\partial G(x,y;x',y')}{\partial y'}\right|_{y'=0}\,dx'\\\\
&=\frac{1}{\pi}\lim_{y\to 0^+}\int_0^\infty g(x')\left(\frac{y}{(x-x')^2+y^2}\right)\,dx'\\\\
&\sim \frac{1}{\pi} g(x)\lim_{y\to 0^+}\int_{-\delta}^\delta \frac{y}{x^2+y^2}\,dx\\\\
&=\frac1\pi g(x) \lim_{y\to 0^+}\left( \left.\arctan(x/y)\right|_{-\delta}^{\delta}\right)\\\\
&=g(x)
\end{align}$$
as expected.
Lastly, we examine $\lim_{x\to 0^+}u(x,y)$. Again, using heuristic formal analysis, we have for "small" $\delta >0$
$$\begin{align}
\lim_{x\to 0^+} \frac{\partial u(x,y)}{\partial x}&=-\lim_{x\to 0^+}\int_0^\infty h(y')\frac{\partial G(x,y;0,y')}{\partial x}\,dy'\\\\
&=\frac{1}{\pi}\lim_{x\to 0^+}\int_0^\infty h(y')\left(\frac{x}{x^2+(y-y')^2}\right)\,dy'\\\\
&\sim \frac{1}{\pi} h(y)\lim_{x\to 0^+} \int_{-\delta}^\delta \frac{y}{x^2+y^2}\,dx\\\\
&=h(y)
\end{align}$$
as expected!
Best Answer
The correct extension is $$\nabla^2G = \delta(\mathbf{x}-\xi) + \delta(\mathbf{x}-\xi_1) - \delta(\mathbf{x}-\xi_2) - \delta(\mathbf{x}-\xi_3)$$ To see why, interpret $G$ as electric potential. Then, the delta functions correspond to point charges, and we have the conditions that (i) the potential is symmetric across $x=0$, and (ii) the potential is zero along $y=0$. The first condition implies that the sign on the delta function at $\xi_1$ should be the same as the sign on the delta function at $\xi$, while the second condition implies that the sign on delta function at $\xi_2$ should be opposite the sign at $\xi$. Then, since $\xi_3$ is the reflection of $\xi_2$ over $x=0$, the delta function there should have the same sign as $\xi_2$ (alternatively, since it is the reflection of $\xi_1$ over $y=0$, it should have the opposite sign from the delta function at $\xi_1$).
Solving for $G$ yields $$G = \ln\frac{1}{|\mathbf{x} - \xi|} + \ln\frac{1}{|\mathbf{x} - \xi_1|} - \ln\frac{1}{|\mathbf{x} - \xi_2|} - \ln\frac{1}{|\mathbf{x} - \xi_3|}$$ and we can verify that $G$ satisfies the required boundary conditions.