You can expand in terms of solutions in one variable if you have a well-posed, fully self-adjoint problem in that direction. For example, suppose $p(x)$ is a real, bounded, continuously differentiable, and strictly positive function on $[a,b]$, and suppose that $q(x)$ is a bounded function on $[a,b]$. Then consider the differential operator $L_{\alpha,\beta}$ defined by
$$
L_{\alpha,\beta}f=-\frac{d}{dx}\left(p(x)\frac{df}{dx}\right)+q(x)f(x)
$$
on the domain $\mathcal{D}(L_{\alpha,\beta})$ consisting of all twice continuously differentiable functions $f$ on $[a,b]$ which satisfy the following endpoint conditions:
$$
\cos\alpha f(a)+\sin\alpha f'(a)=0 \\
\cos\beta f(b)+\sin\beta f'(b)=0.
$$
Then $L_{\alpha,\beta}$ is essentially self-adjoint, meaning that its closure is self-adjoint. It automatically follows that there are eigenvalues $\lambda_1 < \lambda_2 < \lambda_3 < \cdots$ and corresponding normalized eigenfunctions $f_{n}$ such that $L_{\alpha,\beta}f_n=\lambda_n f_n$, and such that $\{ f_{n} \}_{n=1}^{\infty}$ is a complete orthonormal basis of $L^2[a,b]$. (Here $\alpha,\beta$ are used only as a convenient parameterization of the endpoint conditions.)
The inverse operator is compact because of how it takes $L^2$ functions to twice absolutely continuous functions on $[a,b]$. And that's how you end up with a discrete set of eigenvalues $\lambda_1 < \lambda_2 < \lambda_3 < \cdots$ and a set of eigenfunctions $\{ f_{n} \}_{n=1}^{\infty}$ corresponding to these eigenvalues. The eigenfunction spaces are one-dimensional, and all eigenfunctions taken together form a complete orthogonal basis of $L^2[a,b]$. So, every $f \in L^2[a,b]$ may be written as
$$
f(x) = \sum_{n=1}^{\infty}\frac{\langle f,\varphi_n\rangle}{\langle \varphi_n,\varphi_n\rangle}\varphi_n(x) .
$$
This expansion is guaranteed to converge in $L^2[a,b]$, which can be proved using $L^2$ theory only. But, in fact, you'll also have pointwise almost every convergence as well.
The Fourier expansions associated with a Sturm-Liouville operator may switch from being discrete to being "continuous" Fourier integral types of expansions, or even mixed continuous and discrete, when either the interval of consideration is infinite, or when $p$ vanishes at one or both endpoints of the interval $(a,b)$. For example, the radial equation for the hydrogen isotope has both discrete eigenvalues, and has a continuum of spectrum above a certain escape energy. These correspond to bound and unbound states of the atom. The unbound states can have any energy in a continuum, but the bound states must have energies levels that are in a set of discrete values that correspond to specific orbits.
ADDED TO ADDRESS YOUR CASE: To explain more, suppose that you want to solve the Laplace equation
$$
\Delta f = 0, \\
a \le x \le b, \\
c \le y \le d,
$$
subject to boundary conditions on the four edges:
$$
f(x,c) = f_1(x),\;\; f(x,d) = f_2(x),\;\; a \le x \le b, \\
f(a,y) = f_3(y),\;\; f(b,y) = f_4(y),\;\; c \le y \le d.
$$
Then you can solve 4 problems, each one where you specify $f$ to be $0$ on 3 of the 4 edges, and equal to the desired function on the remaining edge. Adding the 4 solutions gives the desired, full solution.
For example, consider the case where $f_2,f_3,f_4$ are $0$, and $f_1$ is not; this is one of the four cases to consider. When you separate variables, the functions in $y$ must vanish at the endpoints $y=c,d$. So you have a nice Sturm-Liouville eigenfunction problem. The standard separation of variables gives
$$
\Delta (XY) = 0 \\
X^{''}Y+XY^{''}=0 \\
\frac{X^{''}}{X}=-\frac{Y^{''}}{Y} \\
\frac{X^{''}}{X}=\lambda,\;\; -\lambda=\frac{Y^{''}}{Y}
$$
The conditions in $x$ do nothing to determine the value of $\lambda$. The values of $\lambda$ are determined by
$$
Y^{''}(y)+\lambda Y(y)=0,\;\; Y(c)=Y(d)=0.
$$
The solutions where $Y(c)=0$ and $Y'(c)=1$ (ignoring $Y(d)$ for the moment) have the form
$$
Y_{\lambda}(y)=\frac{\sin(\sqrt{\lambda}(y-c))}{\sqrt{\lambda}}=\sum_{n=0}^{\infty}(-1)^n\frac{\lambda^n(y-c)^{2n+1}}{(2n+1)!}
$$
NOTE: Adding the normalization $Y_{\lambda}'(c)=1$ forces the above to work in the limiting case as $\lambda\rightarrow 0$. In fact, this general approach will guarantee that the resulting $Y_{\lambda}(y)$ is an entire function of $\lambda$ for every fixed $y\in [c,d]$. $Y_{0}(y)=y-c$ is the correct solution for $\lambda=0$, and it is obtained as a limiting case of the above as $\lambda\rightarrow 0$.
The equation $Y(d)=0$ becomes the eigenvalue equation in $\lambda$:
\begin{align}
&\sin(\sqrt{\lambda}(d-c))=0 \\
&\implies \sqrt{\lambda_n}=\frac{n\pi}{d-c},\;\; n=1,2,3,\cdots, \\
&\implies \lambda_n=\frac{n^2\pi^2}{(d-c)^2},\;\; n=1,2,3,\cdots.
\end{align}
$n=0$ is not a valid case, as was explained above. The eigenfunctions are
$$
\varphi_n(y)=\sin\left(n\pi\frac{y-c}{d-c}\right),\;\; c \le y \le d,\;\; n=1,2,3,\cdots.
$$
$\{ \varphi_n \}_{n=1}^{\infty}$ is a complete orthogonal set of eigenfunctions on $[c,d]$, which is exactly what is needed to solve the problem.
