If we have a rectangular domain $\Omega$ and we are approximating the derivative of $u(x,y)$ by a finite difference $$u_{xy} \approx \frac{(u_{i+1,j+1} + u_{i-1,j-1} - u_{i+1,j-1} - u_{i-1,j+1})}{2dxdy}$$ How would we apply a Neumann boundary condition $\frac{\partial u}{ \partial n} = 0$ on $\partial \Omega$ at the corner of a rectangular domain? E.g. the bottom left green point with u situated at $(0,0)$.
Many thanks.
If a function has $u_x=0$ on $y$-axis, we think of it as symmetric about the $y$-axis (and use this to define its values at the "ghost points" outside of the domain). Same with $u_y=0$ on the $x$-axis. Therefore, if both conditions hold, the function should be treated according to both symmetries, $$u(\pm x, \pm y) = u(x,y)$$ This means all four points marked in green should have the same value of $u$.
As a consequence, $u_{xy}=0$ at the corners, as it is on the boundaries.