Look at a solution $v(x,y)$ of the Neumann problem
\begin{align*}
v_{xx}+v_{yy}& =0\text{ in the unit square} \\
v_{x}(0,y)=v_{x}(1,y)& =0 \\
v_{y}(x,1)& =0 \\
v_{y}(x,0)& =-\pi \cos \pi x.
\end{align*}
Now obtain a solution $u(x,y)$ of your Dirichlet problem as the curve
integral
\begin{equation*}
u(x,y)=\int_{\varphi }\left( -v_{y}\mathrm{\,d}x+v_{x}\,\mathrm{d}y\right) ,
\end{equation*}
where $\varphi $ is any curve in the unit square connecting $(0,0)$ to $(x,y)
$. By the Gauss-Green theorem and the heat equation satisfied by $v$, the
value of this integral is independent of the path. For example,
\begin{equation*}
u(x,y)=-\int_{0}^{x}v_{y}(s,0)\mathrm{\,d}s+\int_{0}^{y}v_{x}(x,t)\,\mathrm{d
}t.
\end{equation*}
From this follows that
\begin{align*}
u_{x}& =-v_{y} \\
u_{y}& =v_{x},
\end{align*}
and hence
\begin{equation*}
u_{xx}+u_{yy}=-v_{yx}+v_{xy}=0.
\end{equation*}
Also the Dirichlet boundary values are easily checked.
Using separation of variables $u = X(x)Y(y)$, we obtain the ODEs
\begin{align}
X'' &= -\lambda X \\
Y'' &= \lambda Y
\end{align}
where we used $-\lambda$ for our separation constant, with associated boundary conditions
\begin{align}
u(L,y) &= 0 \implies X(L) = 0 \\
u_{x}(0,y) &= 0 \implies X'(0) = 0 \\
u_{y}(x,0) - hu(x,0) &= 0 \implies Y'(0) - hY(0) = 0
\end{align}
Solving the ODE in $X$, we find non-trivial solutions only if $\lambda > 0$ which yields
$$X = A \cos \sqrt{\lambda} x + B \sin \sqrt{\lambda} x$$
Now
\begin{align}
X'(0) &= \sqrt{\lambda} B \\
&= 0 \\
\implies B &= 0 \quad \text{(why?)} \\
\therefore X(L) &= A \cos \sqrt{\lambda} L \\
&= 0 \\
\implies \sqrt{\lambda} L &= \frac{(2n + 1) \pi}{2}, \quad n \ge 0 \quad \text{($A \ne 0$ for non-trivial solutions)} \\
\implies \lambda &= \frac{(2n + 1)^{2} \pi^{2}}{4L^{2}}, \quad n \ge 0 \quad(*)
\end{align}
and hence
$$X_{n} = A \cos \left( \frac{(2n + 1) \pi x}{2L} \right)$$
Using the eigenvalue $(*)$ and solving the ODE in $Y$ yields
$$Y = \frac{1}{h} \cdot \frac{(2n + 1) \pi }{2L} \cosh \left(\frac{(2n + 1) \pi y}{2L} \right) + \sinh \left(\frac{(2n + 1) \pi y}{2L} \right)$$
and so the general solution is given by
$$u(x,y) = \sum_{n \ge 0} A_{n} \cos \left( \frac{(2n + 1) \pi x}{2L} \right) \left[ \frac{1}{h} \cdot \frac{(2n + 1) \pi }{2L} \cosh \left(\frac{(2n + 1) \pi y}{2L} \right) + \sinh \left(\frac{(2n + 1) \pi y}{2L} \right)\right]$$
Applying the inhomogeneous condition, we find
\begin{align}
u(x,L) &= u_{0} \\
&= \sum_{n \ge 0} A_{n} \cos \left( \frac{(2n + 1) \pi x}{2L} \right) \left[ \frac{1}{h} \cdot \frac{(2n + 1) \pi }{2L} \cosh \left(\frac{(2n + 1) \pi}{2} \right) + \sinh \left(\frac{(2n + 1) \pi}{2} \right)\right] \\
&= \sum_{n \ge 0} C_{n} \cos \left( \frac{(2n + 1) \pi x}{2L} \right)
\end{align}
where
$$C_{n} = A_{n} \left[ \frac{1}{h} \cdot \frac{(2n + 1) \pi }{2L} \cosh \left(\frac{(2n + 1) \pi}{2} \right) + \sinh \left(\frac{(2n + 1) \pi}{2} \right)\right]$$
You can now solve for the coefficients $C_{n}$ using orthogonality relations. Note that when doing the integrals, a change of variable might help.
Best Answer
The general solution satisfying the boundary conditions $u(x,0)=u(x,1)=0$ is
$$u(x,y)=\sum_{n=1}^{\infty}[A_n\sinh(n\pi x)+B_n\cosh(n\pi x)] \sin(n\pi y),$$
so
$$\partial_x u(x,y)=\sum_{n=1}^{\infty}[n\pi A_n\cosh(n\pi x)+n\pi B_n\sinh(n\pi x)] \sin(n\pi y).$$
Applying the boundary condition $\partial_x u(0,y)=\sin(\pi y)$, the above equation yields
$$\sum_{n=1}^{\infty}n\pi A_n \sin(n\pi y)=\sin(\pi y),$$
which implies $A_1=\frac{1}{\pi}$ and $A_n=0$ for $n\neq 1$. Applying the boundary condition $\partial_x u(1,y)=\sin(3\pi y)$ we obtain
$$[\cosh(\pi)+\pi B_1\sinh(\pi)]\sin(\pi y)+\sum_{n=2}^{\infty}n\pi B_n\sinh(n\pi) \sin(n\pi y)=\sin(3\pi y),$$
which implies $B_1=-\frac{1}{\pi}\coth(\pi),\, B_3=\frac{1}{3\pi}\operatorname{csch}(3\pi) $, and $B_n=0$ for $n\notin\{1,3\}$. Inserting all these results in the expression for $u(x,y)$, we finally obtain
$$u(x,y)=\frac{1}{\pi}[\sinh(\pi x)-\coth(\pi)\cosh(\pi x)]\sin(\pi y)+\frac{1}{3\pi} \operatorname{csch}(3\pi)\cosh(3\pi x)\sin(3\pi y). $$