Suppose you have a differential operator $P(x,y) = P_1(x) + P_2(y)$, and you have a differential equation about $F(x,y)$ $$ P(x,y) F(x,y) = 0 $$ Then one can try to find a basis $\{ F_n(x,y)\}$ for the solution space, where each $F_n$ can be written as (omitting the subscript $n$) $$ F(x,y) = A(x) B(y)$$ Then the equation become $$ B(y) P_1(x) A(x) + A(x) P_2(y) B(y) = 0 $$ Divide both sides by $A(x) B(y)$, we get $$ \frac{1}{A(x)} P_1(x) A(x) + \frac{1}{B(y)} P_2(y) B(y) = 0. $$ The first term is something about $x$, the second term is something about $y$, the only possibility that they cancel is that both terms are constants, hence we try to solve for $$ \frac{1}{A(x)} P_1(x) A(x) = \lambda, \quad \frac{1}{B(y)} P_2(y) B(y) = -\lambda $$ for some $\lambda \in \R$.

Naive Approach If for some $\lambda$, we can solve to get $A_{\lambda,i}(x)$ for $i=1,\cdots, n$, and $B_{\lambda, j}(y)$ for $j =1, \cdots, m$, then we get a collection of functions $$F_{\lambda, i, j} (x, y) = A_{\lambda,i}(x) B_{\lambda, j}(y) $$ and we may hope to express any solution as $$ F(x,y) = \sum_{\lambda} \sum_{i} \sum_{j} C_{\lambda, i, j } A_{\lambda,i}(x) B_{\lambda, j}(y) $$

Problem with the above approach How do I know that these functions $ \{ F_{\lambda, i, j} (x, y) \}$ forms a basis of the solution space?

Consider the equation $$ (\d_x^2 + \d_y^2) F(x,y) = 0 $$ Or the equation $$ (\d_x^4 + \d_y^4) F(x,y) = 0 $$

Consider the example of steady state heat equation in a rectangle. $$ \tag{1} (\d_x^2 + \d_y^2)T(x,y) = 0, \quad 0< x < 1, 0 < y < 1 $$ with boundary condition that $$ \tag{2} T(0, y) = 0, \quad T(1,y) = 0, \quad T(x, 0) = 0, $$ and in addition we impose $$ \tag{3} T(x, 1) = f(x) $$ for some smooth function $f(x)$ on $[0,1]$, such that $f(0)= f(1)=0$.

How to solve this problem? We hope to write the solution in the following form $$ T(x,y) = \sum_{n} C_n X_n(x) Y_n(y) $$ where $X_n(x) Y_n(y)$ forms a basis of the above problem with constraint (1) and (2). Then (3) can be used to fix the coefficients $C_n$.

The separation of variable gives $$ \d_x^2 X(x) = \lambda X(x), \quad X(0) = X(1) = 0 $$ $$ \d_y^2 Y(y) = -\lambda Y(y), \quad Y(0) = 0 $$ What $\lambda$ should we choose? The equation about $X(x)$ indicate $\lambda$ should be negative, we should have $$\lambda = -(n \pi)^2, \quad X_{\lambda}(x) = \sin(n \pi x) $$ And the equation for $Y$ can be solved as $$ Y(y) = \sinh( n \pi y) $$

The general solution for (1) and (2) is $$ F(x,y) = \sum_{n=1}^\infty C_n \sin(n \pi x) \sinh(n \pi y) $$

Finally, we use the boundary condition $f(x) = F(x, 1)$ to get $C_n$ $$ f(x) = \sum_{n=1}^\infty C_n \sin(n \pi x) \sinh(n \pi) $$ Multiply both sides by $\sin(n \pi x)$ and integrate over $x$, we get $$ \int_0^1 f(x) \sin(n \pi x) dx = C_n \frac{1}{2} \sinh(n \pi) $$ hence $$ C_n = \frac{2}{\sinh(n \pi)} \int_0^1 f(x) \sin(n \pi x) dx. $$

I feel a bit lucky that this can be solved.