Can we solve the general eigenfunction of $L_x$?

Consider the orbital angular momentum in QM, labeled by $L$ ($mathbf{L}=mathbf{r}timesmathbf{p}$). In spherical coordinate, the operator can be expressed as:
begin{equation*}
left{begin{aligned}
L_x&=frac{hbar}{mathrm{i}}left(-sinphifrac{partial}{partial theta}-cosphicotthetafrac{partial}{partial phi}right)\
L_y&=frac{hbar}{mathrm{i}}left(cosphifrac{partial}{partial theta}-sinphicotthetafrac{partial}{partial phi}right)\
L_z&=frac{hbar}{mathrm{i}}frac{partial}{partial phi}
end{aligned}right.
,L^{2}=-hbar^{2}left[frac{1}{sin theta} frac{partial}{partial theta}left(sin theta frac{partial}{partial theta}right)+frac{1}{sin ^{2} theta} frac{partial^{2}}{partial phi^{2}}right]
end{equation*}

We know that in general, $[L_z,L_x]=mathrm{i}hbar L_yneq0$, so they do not share the same basis. By solving the eigenquation:
begin{equation*}
left{begin{aligned}
L^2psi&=hbar^2l(l+1)psi\
L_zpsi&=hbar mpsi
end{aligned}right.
end{equation*}
We find the common eigenfunction of $L^2$ and $L_z$ is spherical harmonic funcion $Y_l^m$. But what about $L_x$ and $L_y$? Can we get the general eigenfunction of $L_x$ for state $|l; mrangle$(which means $L_xpsi=hbar mpsi$) using the same method?

I know for certain cases(or each case), the eigenfunction of $L_x$ can be expressed as the linear combination of $Y_l^m$, just use ladder operator to expand the operator in $Y_l^m$ basis. For example, suppose $l=1$, we have $Y_1^{-1},Y_1^0,Y_1^1$ as the basis, so we set:
begin{equation*}
begin{aligned}
Y_1^1=
begin{pmatrix}
1\
0\
0
end{pmatrix}
,
Y_1^0=
begin{pmatrix}
0\
1\
0
end{pmatrix}
Y_1^{-1}=
begin{pmatrix}
0\
0\
1
end{pmatrix}
end{aligned}
end{equation*}
Then, same as what we did for spin, we find the matrix for $L_x$ looks like:
begin{equation*}
begin{aligned}
L_{x}=frac{1}{2}left(L_{+}+L_{-}right)=frac{sqrt{2} hbar}{2}
begin{pmatrix}
0 & 1 & 0 \
1 & 0 & 1 \
0 & 1 & 0
end{pmatrix}
end{aligned}
end{equation*}(Here I omit the calculation), and the eigenstates for $L_x$ are:
begin{equation*}
begin{aligned}
varphi_{hbar}=frac{1}{2}left(begin{array}{c}
1 \
sqrt{2} \
1
end{array}right) varphi_{0}=frac{sqrt{2}}{2}left(begin{array}{c}
1 \
0 \
-1
end{array}right) varphi_{-hbar}=frac{1}{2}left(begin{array}{c}
1 \
-sqrt{2} \
1
end{array}right)
end{aligned}
end{equation*}
Where $varphi_{hbar}$ means the eigenstate with eigenvalue $hbar$. So in this case, the eigenstate can be expressed as the linear combination of $Y_1^m$, namely:
begin{equation*}
begin{aligned}
varphi_{hbar}=frac{1}{2}left(1; sqrt{2}; 1right)
begin{pmatrix}
Y_1^1\
Y_1^0\
Y_1^{-1}
end{pmatrix}
=frac{1}{2}left(Y_1^1+sqrt{2}Y_1^0+Y_1^{-1}right)
end{aligned}
end{equation*}
But we cannot do it in general cases. For example, if the particle is in state $phi=Y_1^0+Y_2^1+Y_4^2$(without normalization), if we want to measure $L_x$, what’s the probability of each value we will get?

I tried to solve the equation like what we did for $Y_l^m$, but I failed. Consider:
begin{equation*}
begin{aligned}
L_x f_l^m(theta,phi)=hbar mf_l^m(theta,phi)Rightarrow frac{hbar}{mathrm{i}}left(-sinphifrac{partial}{partial theta}-cosphicotthetafrac{partial}{partial phi}right)f_l^m(theta,phi)=hbar mf_l^m
end{aligned}
end{equation*}

For PDE, the only way I know is to separate variables: set $f_l^m(theta,phi)=Theta(theta)Phi(phi)$, after plugging in, I found it cannot be solved like usual:
begin{equation}
-frac{1}{Theta}tanthetafrac{mathrm{d}Theta}{mathrm{d}theta}-frac{1}{Phi}cotphifrac{mathrm{d}Phi}{mathrm{d}phi}=mathrm{i}mtanthetafrac{1}{sinphi}
end{equation}
It is not a constant on the right hand side, and I cannot separate it into the sum of two functions, so I don’t know what to do next.

I’v also tried to diagonalize the matrix of $L_x$ in the basis of $L_z$ directly. I find that the matrix of $L_x$ generally looks like:
begin{equation*}
L_{x}=frac{hbar}{2}left(begin{array}{ccccccc}
0 & b_{s} & 0 & 0 & cdots & 0 & 0 \
b_{s} & 0 & b_{s-1} & 0 & cdots & 0 & 0 \
0 & b_{s-1} & 0 & b_{s-2} & cdots & 0 & 0 \
0 & 0 & b_{s-2} & 0 & cdots & 0 & 0 \
vdots & vdots & vdots & vdots & cdots & vdots & vdots \
0 & 0 & 0 & 0 & cdots & 0 & b_{-s+1} \
0 & 0 & 0 & 0 & cdots & b_{-s+1} & 0
end{array}right)
end{equation*}

Where $b_{j} equiv sqrt{(s+j)(s+1-j)}$. Which means the matrix element can be expressed as:$(L_x)_{jk}=frac{hbar}{2}b_{k+1}delta_{j;k+1}+frac{hbar}{2}b_kdelta_{j; k-1}$. Then I write the eigenequation: $L_x|l; mrangle=hbar m|l; mrangle$ as:
begin{equation*}
begin{aligned}
&sum_{k}left(frac{hbar}{2}b_{k+1}delta_{j;k+1}+frac{hbar}{2}b_kdelta_{j; k-1}right)|l; mrangle_k=hbar m|l; mrangle_j\
&Rightarrow frac{1}{2}b_j|l; mrangle_{j-1}+frac{1}{2}b_{j+1}|l; mrangle_{j+1}=m|l; mrangle_j\
&Rightarrow |l; mrangle_{j+1}=frac{2m}{b_{j+1}}|l; mrangle_j-frac{b_j}{b_{j+1}}|l; mrangle_{j-1}
end{aligned}
end{equation*}
If $j>2m+1$ or $j<0$, then $|l; mrangle_j=0$. For the initial condition, we can just set $|l; mrangle_1=1$, after computing every element in the vector, we can then normalize it.

This is a recurrence relation about element of eigenstate of $L_x$, but I cannot solve it. But at least I can comput it.

So my question:

• Can the PDE above be solved to get the general solution of the eigenfunction of $L_x$?(Something looks like $Y_l^m$, maybe some ugly special functions)
• Is there other ways to figure out the general eigenfunction of $L_x$? Or for any $|l; mrangle$, can I calculate the coefficient of $Y_l^m$ quickly? Or solve the series by the recurrence relation I found before?
• Which book can I refer to? I have searched this question on Google, but I found nothing.

Thx for reading my long question. I’m a undergraduate physics student, and I have only learnt QM for a few weeks by myself, so there may be mistakes in the question or misunderstandings for QM. Please point them out if you find, thx!