Introduce $u=X(x)Y(y)$ then applying separation of variables yields
$$\underbrace{\frac{X''}{X}}_{-\lambda_1}+\underbrace{\frac{Y''}{Y}}_{-\lambda_2}=-\lambda$$
Now we have two ODEs, namely $X''+\lambda_1X=0$ and $Y''+\lambda_2Y=0$ both with Neumann B.C. at $0,a$ and $0,b$ respectively. We know the solution to these problems:
$$\begin{cases}
\lambda_{1,0}=0& X_0=\frac{1}{2}\\
\lambda_{1,n}=\frac{n^2\pi^2}{a^2}& X_{n}=\cos(n\pi x/a)\\
\lambda_{2,0}=0& Y_0=\frac{1}{2}\\
\lambda_{2,n}=\frac{m^2\pi^2}{b^2}& Y_{m}=\cos(m\pi x/b)\\
\end{cases}$$
Where $n,m\geq 1$. So the eigenfunctions are
$$u_{n,m}(x,y)=\cos(n\pi x/a)\cos(m\pi x/b)$$
With eigenvalues $\lambda_{m,n}=\lambda_{1,n}+\lambda_{2,n}=\pi^2(n^2/a^2+m^2/b^2)$. We can absorb the case where $m=0$ or $n=0$ by allowing us to plug $m,n=0$ in this equation.
