I think your theorem is a thoroughgoing misconception of the p.d.e. factorization, of, e.g., a spherically symmetric system. Nondimensionalizing all silly constants by absorbing them in the relevant units, you have something like $-\Delta +V(r) -E=0$.
Its eigenvectors are not the product of the eigenvectors of all symmetry generators (operators commuting with the hamiltonian), here, among others, the three $\vec L$s. Instead, recall how this equation's variables separate in plain p.d.e. theory (cf. separation of variables ):
$$
0=-\Delta +V(r) -E\\ = -\frac{1}{r^2}\partial_r ~r^2\partial_r + \frac{1}{r^2} L^2(\theta,\phi) +V(r) -E ~~~\leadsto \\
L^2 (\theta,\phi) = \partial_r ~r^2\partial_r -r^2 V(r) +r^2 E .
$$
This is plain separation of variables: each side of the equation involves different variables, so its eigenvector structure is disjoint. The eigenvectors of the l.h.s. are the spherical harmonics, $Y_{lm}(\theta, \phi)$, with eigenvalues $l(l+1)$ and the eigenvectors of the r.h.side must be
functions of just r, but with the same eigenvalues, i.e.
$$
-\frac{1}{r^2}\partial_r ~r^2\partial_r + \frac{l(l+1) }{r^2} +V(r) =E,
$$
now written in a more familiar form, and with E to be determined, the new eigenvalue, for eigenfunctions $R_{nl}(r)$.
The radial eigenfunctions are inert under $L(\theta,\phi)^2$, but still contain its eigenvalues; and the spherical harmonics do not depend on r, but, of course, are not eigenfunctions of $L_x$, $L_y$, only of $L^2$ and $L_z$, which commute with the hamiltonian as well, so they are good symmetry charges for it.
Now, in other coordinate systems, and for special potentials, like the Hydrogen's you may be more efficient (cf. Pauli's original SO(4) solution of the problem; could do worse than studying this one.), but factorization of p.d.e.s is usually guided by symmetry, as you saw above. You should best consider selected eigenfunctions of the symmetry operators, and utilize those that entangle least with the rest of them.
Finally, in the trivial case where the symmetries commute among themselves, then, of course, the Hilbert space itself factorizes into a tensor product; whose tensor factors, and thus wave function, factor, and are exclusively operated upon by the corresponding eigen-operator, oblivious of the other tensor factors corresponding to the other operators. If this trivial case is the case your instructor discussed, it's hardly salutary to formalize it so impossibly abstractly.
- Response to comment on suppositious example (1).
Let's write it in nondimensionalized units,
$$
0=-E+\vec p^2 +2k \vec p\cdot \vec S .
$$
Without loss of generality, for the purposes of illustrating the problem, take $\vec S= \vec \sigma /2$. Since the problem is manifestly spherically symmetric, we can always rotate the spin to the 3rd (z) direction, without affecting its eigenvalues!
$$
0=-E+\vec p^2 +k p_z \sigma_3~~.
$$
The problem has separated to three decoupled pieces,
$$
0=-E+p_x^2+ p_y^2+( p_z^2 +kp_z\sigma_3 ) .
$$
The first two pieces are scalar, but the third is a 2×2 matrix, so, acting on the space of 2-spinors.
The eigenvalues $E=\vec p^2 \pm kp_z$ are the sum of the eigenvalues of each piece on the right, for the eigenfunctions constant(x), constant'(y) and constant times
$$
\begin{bmatrix} 1\\0\end{bmatrix} ; \qquad \begin{bmatrix} 0\\1\end{bmatrix},
$$
respectively. Note the last eigenfunction hinges on the eigenspinors of $S_3$ only: you couldn't possibly diagonalize all three $S_x,S_y,S_z$ simultaneously. Actually, the first two have completely dropped out of the problem.
In momentum space, which you converted to a meretricious laconism by Fourier-transforming, there is hardly any insight into your suppositious theorem, and you properly invited me to not focus on it. Your hamiltonian is a 2×2 spin matrix and obviously its eigenvectors are 2-spinors.