Solving Simple PDE by Green's Function, Very Confused By Some Mistake

[Metronome]

Suppose I want to solve $u_{xy} = xy$ via Green's Function. This will correspond to the associated PDE $G_{xy} = \delta(x - x_G,\ y - y_G)$, and I want my boundary conditions for this Green problem to be $G(0,\ y) = 0,\ G_x(x,\ 0) = 0$.

The first step logically would actually be to find the adjoint of the differential operator and boundary conditions for $u$ to arrive at the boundary value problem for $G$, but I'll just assume everything is self-adjoint for now, and return to that step at the end. I solve $G_{xy} = \delta(x - x_G,\ y - y_G),\ G(0,\ y) = 0,\ G_x(x,\ 0) = 0$ using a double Laplace Transform. My double-transformed equation is $\hat{\hat G} = \frac{1}{rs}e^{-y_G s - x_G r}$, and unwinding the transforms, my Green's Function turns out to be $G = H[x - x_G]H[y - y_G]$, where $H[\ ]$ is the Heaviside Function.

The solution to the original PDE thus should be $u = \int_0^L \int_0^W x_G y_G H[x - x_G]H[y - y_G]\ dy_G\ dx_G$. Making use of the identity $H[a - b] = 1 - H[b - a]$, I can eliminate the Heaviside Functions by modifying the integration limits, and then straightforward calc computations yield $u = \frac{x^2 y^2}{4}$ (the domain dimensions $L$ and $W$ cancel).

If the operator/boundary conditions are indeed self-adjoint, and if I did everything correctly, then $u_{xy} = xy,\ u(0,\ y) = 0,\ u_x(x,\ 0) = 0$ should yield $u = \frac{x^2 y^2}{4}$. Solving this problem by simple integration and application of the boundary conditions confirms this answer, which in turn strongly confirms the self-adjointness hypothesis as well as my computations so far.

Returning to the beginning, all that remains is to find the adjoint to verify this. I have $u_{xy} = xy$, and I take the inner product of both sides with $G$. So I have to use integration by parts on $\int_0^L \int_0^W u_{xy}G\ dy\ dx$ with respect to $y$. This gives $\int_0^L (u_x G |_{y = 0}^{y = W} - \int_0^W u_x G_y\ dy) dx$. After distributing the outer integral, swapping the order of integration in the double integral, and applying integration by parts again, this time with respect to $x$, I get $\int_0^L u_x G |_{y = 0}^{y = W}\ dx - \int_0^W u G_y |_{x = 0}^{x = L}\ dy + \int_0^L \int_0^W uG_{xy}\ dy\ dx$.

The final term is indeed the original inner product with the differential operator moved from $u$ to $G$. I just have to check the boundary conditions now; I plug $u(0,\ y) = 0$ and $u_x(x,\ 0) = 0$ into $\int_0^L u_x G |_{y = 0}^{y = W}\ dx - \int_0^W u G_y |_{x = 0}^{x = L}\ dy$, and am left with $\int_0^L u(x,\ W)G(x,\ W)\ dx - \int_0^W u(L,\ y)G(L,\ y)\ dy$.

However, my understanding is that I now have to pick two boundary conditions in $G$ that make this expression $0$, and all I have left to use are $G(x,\ W)$ and $G(L,\ y)$. Setting these to $0$ won't give me the $G(0,\ y) = 0$ and $G_x(x,\ 0) = 0$ boundary conditions I expected. This seems to suggest that the operator $\frac{\partial}{\partial y \partial x}$ is actually not self-adjoint. But this doesn't make sense either, because when I assumed self-adjointness, I got the correct answer, and when I use simple integration to solve $u_{xy} = xy$ and apply the boundary conditions in $u$ which would produce the needed boundary conditions in $G$ (they are $u(x,\ W) = 0$ and $u_y(L,\ y) = 0$), I no longer get the solution $u = \frac{x^2 y^2}{4}$ to match my Green's Function solution.

What am I doing wrong here? Is there some simple mistake in my computations, or do I understand the procedure of finding the adjoint operator/boundary conditions incorrectly?