Solving two coupled PDE from a quantum problem

1. Dec 1, 2013

Physteo

Hello Guys.
I have to solve two coupled PDE coming from a quantum physical problem, which possess a cylindrical symmetry. They look like this
$\left\{-\frac{i\hbar}{2M}\left( D^2_\rho + \frac{1}{\rho}D_\rho + D^2_z - \frac{L_z^2}{\hbar^2 \rho^2} \right) + \frac{M\omega^2_\rho \rho^2}{2} + \frac{M\omega^2_z z^2}{2} + g\psi^2_0 - \mu \right\} f^+_\nu = \epsilon_\nu f^-_\nu$
$\left\{-\frac{i\hbar}{2M}\left( D^2_\rho + \frac{1}{\rho}D_\rho + D^2_z - \frac{L_z^2}{\hbar^2 \rho^2} \right) + \frac{M\omega^2_\rho \rho^2}{2} + \frac{M\omega^2_z z^2}{2} + 3g\psi^2_0 - \mu \right\} f^-_\nu = \epsilon_\nu f^+_\nu$

where
$f^{+,-}_\nu = f^{+,-}_\nu(\rho, z, \phi)$

$D_x \equiv \frac{\partial}{\partial x}$

$\psi_0^2 = \frac{\mu}{g} \left(1 - \frac{\rho^2}{R^2_\rho} - \frac{z^2}{R^2_z}\right)$

$L_z$ is the projection of the angular momentum operator, which operates like this on his eigenfunctions:
$L_z (\frac{\exp{(im\phi)}}{\sqrt(2\pi)}) = \hbar m \frac{\exp{(im\phi)}}{\sqrt(2\pi)}$

and $\hbar, M, \omega_\rho, \omega_z, R_\rho, R_z, \epsilon_\nu$ and $\mu$ are constants.

Actually, I have to find the energies $\epsilon$.
I would like to use the fact that I know how the operator Lz acts, but also, I know how the operator on the z coordinate acts. I would like to end up with only a radial equation on ρ.
To be clearer, I rewrite the equations like this (grouping operators acting on different variables, and plugging the expression for $\psi_0$):

$\left\{-\frac{i\hbar}{2M}\left( D^2_\rho + \frac{1}{\rho}D_\rho \right) + \frac{M\omega^2_\rho \rho^2}{2} -\frac{\rho^2}{R_\rho} \right\} f^+_\nu +\left\{ -\frac{i\hbar}{2M} D^2_z + \frac{M\omega^2_z z^2}{2} -\frac{z^2}{R_z}\right\}f^+_\nu -\frac{i\hbar}{2M} \left\{-\frac{L_z^2}{\hbar^2 \rho^2} \right\}f^+_\nu= \epsilon_\nu f^-_\nu$

And here comes my question. I know how Lz acts on his eigenfunctions, but I also know how the operator on z acts on his eigenfunctions ( because it is the 1D quantum oscillator).
Then, am I allowed to write the solution in the following way?
$f^{+,-}_\nu = \sum_{m} \sum_{n} f^{+,-}_{n,m,\nu}(\rho) \psi_n(z) \frac{\exp{(im\phi)}}{\sqrt(2\pi)}$

where
$\psi_n(z)$ are the eigenfunctions of the harmonic oscillator?

If this is possible, the equation will be simplified in a 1D equation in ρ, replacing the harmonic operator in z, and the operator Lz with their eigenvalues.
Thanks for the attention.