Salmone said:
I didn't mean this, what I mean was: if we are using the method of "separable Hamiltonians" don't we have to prove that the Hamiltonian, not the SE, splits in a sum of terms commuting with each other?
That is too restrictive. The thing to note is that if you have a linear differential operator on the form
$$
\hat L = \hat L_x + f(x) \hat \Lambda
$$
where ##\hat \Lambda## is a Sturm-Liouville operator in some variable ##y## and ##\hat L_x## is a linear differential operator, then you can expand any function ##g(x,y)## in the eigenfunctions ##u_n(y)## of ##\hat\Lambda##:
$$
g(x,y) = \sum_n g_n(x) u_n(y).
$$
For the term ##g_n(x) u_n(y)##, ##\hat L## becomes
$$
\hat L g_n(x)u_n(y) = u_n(y) [\hat L_x + \lambda_n f(x)] g_n(x).
$$
If ##\hat L_x + \lambda_n f(x)## is a Sturm-Liouville operator, then any function ##g_n(x)## can be expanded in its eigenfunctions ##v_{nm}(x)## as
$$
g_n(x) = \sum_m G_{nm} v_{nm}(x).
$$
It is clear that the products ##\psi_{nm}(x,y) = v_{nm}(x) u_n(y)## are eigenfunctions of the full differential operator ##\hat L## and that any function ##g(x,y)## can be expressed as a linear combination
$$
g(x,y) = \sum_{n,m} G_{nm} \psi_{nm}(x,y).
$$
Note that the operator ##\hat L_x + \lambda_n f(x)## is
not generally independent of the eigenfunctions of the ##\hat\Lambda##, but different ##n## lead to different operators in the ##x## direction through different ##\lambda_n##.