The representation of angular momentum (of $L^2$ and $L_z$) in the space is given by
$$ \langle\vec{r}|L_z|\psi\rangle = \frac{\hbar}{i}\frac{\partial}{\partial\phi}\langle\vec{r}|\psi\rangle $$
and $$ \langle\vec{r}|L^2|\psi\rangle = -\hbar^2\left[\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\phi^2} + \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)\right]\langle\vec{r}|\psi\rangle $$
I want to solve the Zeeman's effect problem, consisting of a Hamiltonian of the form
$$ H = \frac{L^2}{2\mu R^2} - \frac{eBL_z}{2\mu}, $$ where $\mu,R,e,B$ are constants. I have to solve the equation
$$ H|\psi\rangle = E|\psi\rangle. $$
Clearly the eigenstates must be eigenstates of $L^2$ and $L_z$ simultaneously. I write this as $|\psi\rangle = |\ell, m \rangle$. Taking into account that we define
$$ L^2|\ell, m\rangle = \hbar^2\ell(\ell+1)|\ell, m\rangle, \quad\quad L_z|\ell, m\rangle = \hbar m |\ell, m\rangle. $$
Then we can find the eigenenergies: $$ E_{\ell, m} = \frac{\hbar^2\ell(\ell+1)}{r\mu R^2}-\frac{e\hbar B}{2\mu}m $$
I want to find the eigenstates, therefore I do the following: $$ \langle\theta, \phi | H |\ell,m\rangle = E_{\ell,m}\langle\theta, \phi |\ell,m\rangle. $$
This differential equation is difficult to solve. Is my approach correct to find the eigenstates?
Thank you for your help.