Consider this problem: Let a 2D box of side length $L$ lie in the first quadrant of the $x$-$y$ plane with one corner at the origin. The vertical walls have arbitrary (but specified) potentials $\phi_1(0,y)$ and $\phi_2(L,y)$, and the horizontal walls have potentials $\phi_3(x,0)$ and $\phi_4(x,L)$. Note that these boundary conditions are functions, not constants. Let $\phi(x,y)$ be defined inside the box and obey $\nabla^2\phi(x,y)=0$. Find $\phi(x,y)$.
Problems like this are often given in physics, but with extremely simple boundary conditions that allow one to set many constants to zero (on the separated solutions) and restrict allowed values of the separation constant. But does the method work for all possible boundary conditions? (Let's consider only $\phi$ boundary functions that have a finite number of discontinuities—nothing pathological.)
In particular, I'm not sure how we will enforce the final two boundary conditions. We can use the sin/cos solutions only to make our solution match one pair of parallel walls, but for the other pair, we are left with exponential functions or sinh/cosh functions, which are not orthogonal. How does one proceed?
For linear equations, which is all one usually considers with separation of variables since otherwise solutions can't be added, there is a rather cheap-looking procedure that one employs to deal with this sort of issue: we exploit the linearity of the equation a bit more to split the boundary conditions into something easier to deal with. Namely, let $u_1$ solve the problem $$ \begin{cases} \nabla^2 u_1 = 0 & 0<x<L,0<y<L \\ u_1(0,y) = \phi_1(0,y) & 0<y<L \\ u_1 = 0 & \text{on the other 3 boundaries} \end{cases}, $$ and similarly $u_2$ solves $$ \begin{cases} \nabla^2 u_2 = 0 & 0<x<L,0<y<L \\ u_2(L,y) = \phi_1(L,y) & 0<y<L \\ u_2 = 0 & \text{on the other 3 boundaries} \end{cases}, $$ and likewise $u_3$ and $u_4$. Then taking $\phi = u_1+u_2+u_3+u_4$ solves the original problem.
Solving the 4 new problems is much easier, of course: we just need a sine series for $\phi_i$, which automatically satisfies the two zero boundary conditions on the two boundaries opposite one another, and a $\frac{\sinh{kx}}{\sinh{kL}}$ or $\frac{\sinh{k(x-L)}}{\sinh{kL}}$ for the other direction to get the zero on the other boundary.