Separation of Variables and Linear PDEs Separation of variables is a powerful method which comes to our help for finding a closed form solution for a linear partial differential equation (PDE). For example, we all know that how the method works for the two dimensional Laplace equation in Cartesian Coordinates
$$\nabla^2 \phi(x,y) = \partial^{2}_{x}\phi(x,y) + \partial^{2}_{y}\phi(x,y) = 0 \tag{1}$$
The steps are
$1$. Consider a solution of the form $\phi(x,y)=X(x)Y(y) \tag{2}$
$2$. Put it into the equation, assuming $X(x) \ne 0,$ and $Y(y) \ne 0$, to get 
$$\dfrac{\text{D}^2 X}{X}+\dfrac{\text{D}^2 Y}{Y}=0 \tag{3}$$
where $\text{D}$ is the differential operator.
$3$. Observe that the above functional equation $(3)$ is possible if and only if
$$\dfrac{\text{D}^2 X}{X} = -\dfrac{\text{D}^2 Y}{Y} = \lambda \tag{4}$$
where $\lambda$ is some constant.
$4$. Conclude that $X$ and $Y$ should satisfy the following ordinary differential equations (ODEs)
$$\begin{align}
\left[ \text{D}^2 - \lambda \right]X &= 0 \\
\left[ \text{D}^2 + \lambda \right]Y &= 0
\end{align} \tag{5}$$
However, I always had some questions about this method which is not addressed in the elementary books. Here are my questions
$1$. What are the restrictions of this method? I mean when it is not going to work! or equivalently When a linear PDE cannot have a separable solution of the form $(2)$?
$2$. Has the method been used for system of linear PDEs? If Yes, would you please give an example.
 A: The point of separation of variables is not just to get some solution, but to get a general solution, which can be used to produce a solution for any initial condition. Consider the case of the initial value problem $u'(t) = u(t)$, $u(0)=1$. You solve this problem by first writing down the general solution $u(t) = C e^t$, and then plugging in the initial condition to find that $C=1$. 
Separation of variables is similar. Consider the heat equation, $u_t = u_{xx}$, on $\{(t,x): 0<t<\infty, 0<x<1\}$, subject to the Dirichlet boundary conditions $u(0,t) = u(1,t) = 0$ for all $t$ and the initial condition $u(x,0) = f(x)$. The steps are as follows: first find all non-zero functions of the form $u(x,t) = X(x)T(t)$ which satisfy the equation $u_t = u_{xx}$ and the boundary conditions $u(x,0) = u(x,1) = 0$. This amounts to demanding that $XT' = X''T$ and that $X(0) = X(1) = 0$. Using the sorts of tricks you mentioned (if you haven't done this example yet, try it!), we can conclude that $X$ must be of the form $X_n(x) = B_n \sin (\pi n x)$, for any $n \in \mathbb N$ and any real number $B_n$, and $T_n = e^{- (\pi n)^2 t}$. We then use the fact that the equation and the boundary conditions are linear, which is sometimes called superposition by fancy folks. That is, if $u_1$ and $u_2$ solve the equation and satisfy the boundary conditions, and $a,b \in \mathbb R$, $a u_1 + b u_2$ is a solution.
It follows that, provided that the $B_n$ decay rapidly enough that the sum converges in some reasonable sense, the function $u(x,t) = \sum\limits_{n=1}^\infty B_n \sin(\pi n x) e^{-(\pi n)^2 t}$ solves the equation and satisfies the boundary condition. Here is where I think an answer to at least your first question comes in. We have that $u(x,0) = \sum\limits_{n=1}^\infty B_n \sin(\pi n x)$. It turns out that every reasonable function (say piecewise continuous) can be written in this form (with the convergence interpreted suitably). You can find the coefficients using techniques you are probably familiar with.
The reason this worked has to do with steps 3 and 4 in your question. After you separate variables, you find that $X$ must solve the boundary value problem $X'' = \lambda X$, $X(0) = X(1) = 0$. This is called a Sturm-Liouville problem, and there is a very well-developed theory for these problems, see Wikipedia. The point is that, under suitable assumptions, you are guaranteed that the solutions to the Sturm Liouville problem form a basis for $L^2([0,1])$. That is, you get a sequence of eigenvalues $\lambda_n$ and solutions of the boundary value problem $X_n$ such that for any (reasonable) function $f$ on $[0,1]$, $f(x) = \sum\limits_{n=1}^\infty B_n X_n(x)$. Here the infinite sum means that 
$$\lim\limits_{n \to \infty} \int\limits_0^1 |\sum\limits_{i=1}^n B_n \sin(\pi n x) - f(x)|^2 dx = 0,$$
but stronger forms of convergence can be shown for certain cases.
In the example case, the $X_n$ were the functions $X_n(x) = \sin (\pi n x)$ and the $\lambda_n$ were the numbers $- (\pi n)^2$. 
So, I think, the answer to your first question is whenever the resulting Sturm-Liouville boundary value problem is regular, and therefore gives a nice set of solutions which can be used to reconstruct the initial condition. 
A: For linear equations, the technique of separation of variables is used to find all separated solutions of the form $X_1(x_1)X_2(x_2)\cdots X_n(x_n)$. You find them all if your equation can be separated.
If you make a change of variables, then you will generally find a different set of separated solutions. For example, you might separate $X(x)Y(y)$, or instead $R(r)\Theta(\theta)$ where $r=\sqrt{x^2+y^2}$, etc. The separated solutions you end up with using a different set of coordinates will be different because the separated solutions $u(x,y)=X(x)Y(y)$ that vanish at some $x$ vanish for all $y$ at that $x$; so the zero sets are coordinate-aligned. The solutions $u(r,\theta)=R(r)\Theta(\theta)$ will vanish on circles and or rays starting at the origin.
Note: just because you can separate in one set of coordinates does not mean you can separate in another.
The Laplacian separates in Cartesian coordinates, and one contributing factor is that there are no mixed derivatives. Coordinate changes that introduce mixed terms most often result in equations that are no longer separable. For example, try $u=X(x)Y(y)$ in the following
$$
           u_{xx}+ 2u_{xy} +u_{yy} = 0 \\
          X''Y + 2 X'Y' + XY'' = 0 \\
          \frac{X''}{X}+2\frac{X'}{X}\frac{Y'}{Y}+\frac{Y''}{Y}=0.
$$
For the Laplacian, avoiding mixed terms essentially requires that you use orthogonal coordinate systems to separate variables, but that's not always enough to be able to separate. The spherical coordinate system and the cylindrical coordinate system are orthogonal because the coordinate surfaces $r=A$, $\theta=B$, $\phi=C$ are mutually orthogonal where they intersect. So, when you separate variables you don't end up with mixed derivative terms. As I recall, there are 24 orthogonal coordinate systems where the Laplacian separates. This limits you, because you can only impose conditions on surfaces where a coordinate is constant, but there are enough systems that significant problems could be solved exactly a century ago, enough to develop Quantum Mechanics, for example.
If you're going to use the technique on a system, then all of the equations need to separate, and maybe in the same coordinate system.
The basic reason you have enough separated solutions to build up a full solution for classical operators is that they are formally selfadjoint. Symmetry is a big part of what makes separation of variables work to give you a full range of solutions. For example,
$$
          \int_{\Omega} \nabla^2 f g dx = \int_{\Omega} f \nabla^2 g dx + \mbox{eval terms}
$$
Fortunately, symmetry and Physics combine in a powerful way; and symmetry and math combine in a powerful way. It's difficult to do much in either field without symmetry or something close to it.
A: *

*There might be solutions which can not be written as $\phi(x,y)=f(x)\, g(y)$, e.g. $\phi(x,y) = f(x+y)$.

*The separation turns each PDE of $n$ variables into a system of $n$ ODEs (one variable each). If common solutions exist, they should assemble to a valid solution of the PDEs as well. How likely this is I do not know. Maybe someone has a nice example?
