Can we solve the general eigenfunction of $L_x$? Consider the orbital angular momentum in QM, labeled by $L$ ($\mathbf{L}=\mathbf{r}\times\mathbf{p}$). In spherical coordinate, the operator can be expressed as:
\begin{equation*}
\left\{\begin{aligned}
  L_x&=\frac{\hbar}{\mathrm{i}}\left(-\sin\phi\frac{\partial}{\partial \theta}-\cos\phi\cot\theta\frac{\partial}{\partial \phi}\right)\\
  L_y&=\frac{\hbar}{\mathrm{i}}\left(\cos\phi\frac{\partial}{\partial \theta}-\sin\phi\cot\theta\frac{\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_y\neq0$, so they do not share the same basis. By solving the eigenquation:
\begin{equation*}
\left\{\begin{aligned}
 L^2\psi&=\hbar^2l(l+1)\psi\\
 L_z\psi&=\hbar m\psi
\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\; m\rangle$(which means $L_x\psi=\hbar m\psi$) 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}\; 1\right)
\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(-\sin\phi\frac{\partial}{\partial \theta}-\cos\phi\cot\theta\frac{\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}\tan\theta\frac{\mathrm{d}\Theta}{\mathrm{d}\theta}-\frac{1}{\Phi}\cot\phi\frac{\mathrm{d}\Phi}{\mathrm{d}\phi}=\mathrm{i}m\tan\theta\frac{1}{\sin\phi}
\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_k\delta_{j\; k-1}$. Then I write the eigenequation: $L_x|l\; m\rangle=\hbar m|l\; m\rangle$ as:
\begin{equation*}
\begin{aligned}
 &\sum_{k}\left(\frac{\hbar}{2}b_{k+1}\delta_{j\;k+1}+\frac{\hbar}{2}b_k\delta_{j\; k-1}\right)|l\; m\rangle_k=\hbar m|l\; m\rangle_j\\
 &\Rightarrow \frac{1}{2}b_j|l\; m\rangle_{j-1}+\frac{1}{2}b_{j+1}|l\; m\rangle_{j+1}=m|l\; m\rangle_j\\
 &\Rightarrow |l\; m\rangle_{j+1}=\frac{2m}{b_{j+1}}|l\; m\rangle_j-\frac{b_j}{b_{j+1}}|l\; m\rangle_{j-1}
\end{aligned}
\end{equation*}
If $j>2m+1$ or $j<0$, then $|l\; m\rangle_j=0$. For the initial condition, we can just set $|l\; m\rangle_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\; m\rangle$, 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!
 A: The problem lies in your understanding of the operator algebra and Hilbert spaces not so much on your calculus I think. Let us break it down:

Can the PDE above be solved to get the general solution of the eigenfunction of Lx?(Something looks like Yml, maybe some ugly special functions)

The solutions are given by linear combinations of the $Y_\ell^m$, which are obtained by expressing the eigenvectors of $L_x$ in terms of those of $L^2$ and $L_z$. These are unique for a set of given initial conditions, so these are THE solution. The equation they satisfy (the one for $L_x$) is not particularly useful so the solutions don't have a name if that is what your are looking for.
Notice that the coordinate system from the start picks $z$ as a special axis, but your axes can be rotated and in reality $L_x, L_y$ and $L_z$ stand on equal footing. If the question is whether the eigenvectors will be simple harmonics, well ... no, because, given a set of coordinates as chosen, that would mean that $L^2$, $L_z$ and $L_x$ commute with each other but they don't, that is why the spherical harmonics since they are already eigenfunctions of $L^2$ and $L_z$ cannot satisfy the equation for $L_x$, but a linear combination of them can, after all they are a basis for the Hilbert space in question.

Is there other ways to figure out the general eigenfunction of Lx? Or for any |lm⟩, can I calculate the coefficient of Yml quickly? Or solve the series by the recurrence relation I found before?

Avoid going to a specific coordinate representation. Namely, use the operator formalism, it is more general and quicker. If you already know how $L_x$ is written by using ladder operators you can write a matrix simply and diagonalize it. This also provides the transformation from the basis of $|\ell\;m_z\rangle$ to the basis $|\ell\;m_x\rangle$. Your series must lead to the same results if it was done correctly, however through a lengthier path.

Which book can I refer to? I have searched this question on Google, but I found nothing.

You seem to be missing the Hilbert space and linear algebra connections so I would recommend you look for math books on those to get a grip on linear operators and bases. Then you can move on to any standard textbook of quantum mechanics such as Sakurai's (both basic and Advanced). I personally also like Cohen-Tannoudji's two volume set on quantum mechanics.
A: First, look at the $L=1$ case:
$$ Y_1^1 = N(-x - iy) $$
$$ Y_1^0 = \sqrt 2 Nz $$
$$ Y_1^{-1} = N(x - iy) $$
where $N = \frac 1 2 \sqrt{\frac{3}{2\pi}}$.
To get the eigenvalues in $x$, you need to do a coordinate change:
$$ (x, y, z) \rightarrow (y, z, x)$$
so in that basis (I'll call it $F$):
$$ F_1^1 = N(-y - iz) $$
$$ F_1^0 = \sqrt 2 Nx $$
$$ F_1^{-1} = N(y - iz) $$
are the special functions that are eigenvalues of $L_x$. (Note that they are neither more or less ugly than the eigenvalues of $L_z$.)
You can easily verify that the $F_1^m$ satisfy the eigenvector relationships with respect to the $Y_1^m$ you have derived.
For all $l$, note that coordinate substitution is just a 120-degree rotation about the unit vector:
$$ \hat n = \frac 1 {\sqrt 3}(\hat x + \hat y +\hat z) $$
which can be implemented with the Wigner-D matrices. They rotate spherical harmonics to new coordinates:
$$ D^j_{m'm}(\alpha, \beta, \gamma) = \langle jm'|R(\alpha, \beta, \gamma)|jm\rangle$$
where the $(\alpha, \beta, \gamma)$ are the Euler angles in the Z-Y-Z convention. (Since I only rotate for Mars landings, I don't use the notoriously unstable Euler angles, so I can't help you).
Note that each $j$ is an irreducible representation of SO(3), so that any rotation of the $2j+1$ spherical harmonics is closed (hence there is no $j'$): that is, the eigenfunctions of $(L^2, L_x)$ with eigenvalue $(l, m')$ can be completely expressed in terms of eigenfunctions of $(L^2, L_z)$ with eigenvalues $l$ and $m \in (-l, \ldots, +l)$.
