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 point $\xi$ should be considered fixed in this context: we plant it somewhere and then look for Green's function, considered as a function of $x$ only. So, no such things as $\xi_x$ are needed in the calculations.
Also, the signs you chose are incorrect. You need the symmetry properties $$G((-a,b),(c,d)) = G((a,b),(c,d))\tag1$$ $$G((a,b),(-c,d)) = -G((a,b),(c,d))\tag2$$ Even reflection gives zero normal derivative, odd reflection gives zero value.
So, the formula is
$$G(\textbf x , \xi)=\frac1{2 \pi} \ln |\textbf x - \xi | + \frac1{2 \pi} \ln |\textbf x - \xi_1 |-\frac1{2 \pi} \ln |\textbf x - \xi_2 | -\frac1{2 \pi} \ln |\textbf x - \xi_3 |$$
which can be expanded in coordinates as
$$ G((a,b),(c,d))=\frac1{ 4\pi}\left( \ln ((a-c)^2+(b-d)^2) + \ln ((a+c)^2+(b-d)^2) \\ - \ln((a-c)^2+(b+d)^2) - \ln ((a+c)^2+(b+d)^2)\right) $$ This function satisfies (1) and (2), and therefore satisfies the required boundary conditions.
and $$G((a,b),(c,d)) = -G((a,b),(c,d))$$