Consider the following non-homogeneous Poisson equation in a rectangular domain $x\in [-L_x/2,L_x/2]$ and $y\in[0,L_y]$, $$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f(x,y)=a\sin(bx),\quad u(x,y=0)=0,\quad \frac{\partial u}{\partial x}\bigg\vert_{x=+ L_x/2}=g(x,y)=c,\quad \frac{\partial u}{\partial x}\bigg\vert_{x=- L_x/2}=g(x,y)=c,\quad \frac{\partial u}{\partial y}\bigg\vert_{y=L_y}=0,$$ where, $a,b,$ and $c$ are real positive constants.
How to approach this problem? On making the boundary conditions homogeneous by using $u=u_0+v$, where $\partial u_0/\partial_x\vert_{x=\pm L_x/2}=c$, the eigenvectors of the Laplacian ($v_{mn}$), with Neumann conditions at $x=\pm L_x/2$ and $y=L_y$, and a Dirichlet condition at $y=0$, may be formally used as a basis to find the solution to the inhomogeneous PDE in $v$,
$$\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}=f(x,y)=a\sin(bx),\quad v(x,y=0)=0,\quad \frac{\partial v}{\partial x}\bigg\vert_{x=+ L_x/2}=0,\quad \frac{\partial v}{\partial x}\bigg\vert_{x=- L_x/2}=0,\quad \frac{\partial v}{\partial y}\bigg\vert_{y=L_y}=0.$$ However, one does end up with the trivial solution for the eigenvectors, $v_{mn}=0$, in that case.
A numerical solution yields a profile, which "looks like" a sinusoid$\sim\sin\big(\frac{\pi y}{2L_y}\big)$ along the $y$ direction. The corresponding system with $f=0$ does have the simple solution $D\sinh\big(\frac{\pi x}{2L_y}\big)\sin\big(\frac{\pi y}{2L_y}\big)$, with $D$ determined by boundary conditions at $x=\pm L_x/2$.
Even if an exact analytical solution is not possible, analytical approximations are equally welcome.