To preface, I am very new to PDEs and not overly familiar with hyperbolic trig so forgive me if there are obvious mistakes - I have tried my best to research this problem before asking this question. Additionally, I am not an applied mathematician so intuition and applications don't really help my understanding. The question: $$ u_{xx} + u_{yy} = 0, \;\; 0<x<1,\;\; 0<y<1 \\ u(x,0) = f(x), \;\; u(x,1) = 0 \\ u_x(0,y) = u_x(1,y) = 0 $$ Where $f(x)$ is an integrable function defined on $0\leq x\leq 1$ that has vanishing derivative at $x=0$ and $x=1$.
My Solution:
Seek a solution of the form $u(x,y) = X(x)Y(y)$, then the PDE becomes:
$$ X''Y + XY'' = 0 \implies \frac{X''}{X}=-\frac{Y''}{Y} = \lambda $$
Implement homogeneous boundary conditions:
$$ u_x(0,y) = X'(0)Y(y) = 0 \implies X'(0) = 0 \\ u_x(1,y) = X'(1)Y(y) = 0 \implies X'(1) = 0 \\ u_x(x,1) = X'(x)Y(1) = 0 \implies Y'(1) = 0 $$
Now, consider the ODE for $X(x)$:
$$ X'' - \lambda X = 0 $$
We require $\lambda < 0$ so introduce $\omega$ such that $\lambda = - \omega^2$. Thus, we get general solution,
$$\begin{align} X &= A\sin(\omega x) + B\cos(\omega x)\\ \implies X' &= A\omega \cos(\omega x) - B\omega \sin(\omega x)\\ \\ X'(0) &= A\omega \implies A = 0, \;\; X'(1) = -B\omega \sin(\omega) = 0 \implies \sin(\omega) = 0 \\ \end{align}$$
Thus, we get the solution for $X(x)$:
$$ X_n(x) = B_n \cos(\omega_n x), \;\; \lambda_n = -\omega_{n}^2, \;\; \omega_n = \pi n,\;\; n = 1, 2, 3,... $$
Now, the ODE for $Y(y)$ becomes:
$$ Y_n'' - \omega_{n}^2 Y_n = 0 $$
Which has general solution,
$$ Y_n = C_n \sinh(\pi n y) + D_n \cosh(\pi n y) $$
And, this is where I am stuck... Implementing the homogeneous boundary condition on $Y_n$, we get
$$ Y_n(1) = C_n \sinh(\pi n) + D_n \cosh(\pi n) = 0 $$
But, it seems to me that the only way to satisfy this is with the trivial solution $C_n = D_n = 0$. So, my question is, am I missing something with how to implement this boundary condition? Or, have I done something wrong earlier in the question?
I first want to note that there is a solution you missed corresponding to $\lambda=0$: $$ X_0(x) = B_0, \ \ \ \ Y_0(y) = C_0y + D_0. $$
I am going to change notation slightly from yours and absorb the constant $B_n$ into $C_n$ and $D_n$. Since they are all arbitrary at this point, this is no issue. The general solution is then $$ \begin{aligned} u(x,y) &= \sum_{n=0}^\infty X_n(x)Y_n(x) \\ &= C_0y + D_0 + \sum_{n=1}^\infty \cos n\pi x\cdot\left(C_n\sinh n\pi y + D_n\cosh n\pi y\right). \end{aligned} $$
This function satisfies the boundary conditions at $x=0,1$ from construction, so we now apply the boundary condition at $y=0,1$ to find the final solution.
$y=1$:
Applying this boundary condition leads to the following relations:
$$ \begin{aligned} C_0 + D_0 &=0 \\ C_n\sinh n\pi + D_n\cosh n\pi &=0, \ \ \ \ n=1,2,\dots \end{aligned} $$
$y=0$:
Applying this boundary condition leads to the following relation:
$$ \begin{aligned} u(x,0) = D_0 + \sum_{n=1}^\infty D_n\cos n\pi x = f(x). \end{aligned} $$
We must now determine $D_n$ from the Fourier cosine series of $f$ on $[0,1]$, but notice that the equation for $u(x,0)$ is already written as a Fourier cosine series, so the $D_n$'s are exactly the coefficients of the Fourier cosine series of $f$ on $[0,1]$. In particular, we have
$$ \begin{aligned} D_0 &= \int_0^1 f(x)dx \\ D_n &= 2\int_0^1f(x)\cos n\pi x \ dx, \ \ \ \ n=1,2,\dots \end{aligned} $$
We have now determined the values of $D_n$ for all $n$, and we substitute them into our previously obtained equations to determine $C_n$.
$$ \begin{aligned} C_0 + D_0 &=0 \\ C_n\sinh n\pi + D_n\cosh n\pi &=0, \ \ \ \ n=1,2,\dots\\ &\implies \\ C_0 &= -D_0 \\ C_n &= -D_n\coth n\pi, \ \ \ \ n=1,2,\dots \end{aligned} $$
We can now write our final solution: $$ \begin{aligned} u(x,y) &= a_0(1-y) + \sum_{n=1}^\infty a_n\cos n\pi x\cdot\left(\cosh n\pi y - \coth n\pi \cdot\sinh n\pi y\right), \ \ \ \ \text{where} \\ a_0 &= \int_0^1 f(x)dx, \\ a_n &= 2\int_0^1 f(x)\cos n\pi x \ dx, \ \ \ \ n=1,2,\dots \end{aligned} $$