有大佬知道这个方程怎么解吗

这是个数学物理问题，解就是球谐函数，具体步骤可以找一本数学物理方法的教材。
Mathematica有用来表示球谐函数的特殊函数。

这个问题并不容易，因为截止13.0版本，DEigensystem还给不了通解。（算几个具体解还是可以的。）适当加入人工分析或许可以，不过我现在有事要离开一会儿了。

事实上如果你仔细读SphericalHarmonicY的文档就会发现文档编写者其实回避了相关的讨论。Mathematica如果可以漂亮地把它解出来，“属性与关系”一节不可能一个字不提。

来说说怎么用Mathematica解这个方程吧。如前所述，目前（截止13.0）Mathematica还没法一行代码解出这个特征值问题的通解，我们得引入少量的人工分析。首先，方程中对φ的导数显然可以使用有限域傅立叶变换消掉：
(* 我知道教科书在讨论此问题时一般会用分离变量法，但我个人因为这里（http://www.zhihu.com/question/412203550/answer/1388426673）所说的理由，并不喜欢分离变量 *)
eq = With[{Y = Y[θ, φ]},
Laplacian[Y, {r, θ, φ}, "Spherical"] + l (l + 1) Y == 0] /. r -> 1 //
Simplify
(* finiteFourierTransform的定义可以在这帖 mathematica.stackexchange.com/q/155817/1871 找到 *)
teq = finiteFourierTransform[eq, {φ, -Pi, Pi},
m] /. {Y[θ, π] -> Y[θ, -Pi],
Derivative[0, 1][Y][θ, π] -> Derivative[0, 1][Y][θ, -Pi]} /.
HoldPattern@finiteFourierTransform[a_, __] :> a /. Y -> (y@# &)
(*
l (1 + l) y[θ] - m^2 Csc[θ]^2 y[θ] + Cot[θ] y'[θ] + y''[θ] == 0
*)
(* 变换所得的teq是个常微分方程，可以使用DSolve求解： *)
tsol = DSolveValue[teq, y[θ], θ]
(* C[1] LegendreP[l, m, Cos[θ]] + C[2] LegendreQ[l, m, Cos[θ]] *)
(* 但这个结果相较于球谐函数SphericalHarmonicY的定义，多了一个LegendreQ，要怎么样，或者说，要以什么理由把这项给拿掉呢？别忘了，这个问题还隐含着“解在整个定义域内有界”这一限制，而在x=±1时LegendreQ[l,m,x]是无界的 *)
(* 这个结果也可以用Mathematica辅助计算导出。只要使用12.2引入的FunctionDiscontinuities：*)
FunctionDiscontinuities[LegendreQ[l, m, Cos[θ]], θ] //
FullSimplify // LogicalExpand
(* (l >= m && Cos[θ] >= 1) || (l >= m &&
Cos[θ] <= -1) || (Cos[θ] >= 1 && l + m < 0) || (Cos[θ] >= 1 &&
l ∉ Integers) || (Cos[θ] >= 1 &&
m/2 ∉ Integers) || (l + m < 0 &&
Cos[θ] <= -1) || (Cos[θ] <= -1 &&
l ∉ Integers) || (Cos[θ] <= -1 && m/2 ∉ Integers) ||
Cos[θ] == -1 || Cos[θ] == 1 *)
(* 级数展开也是可以的： *)
Series[tsol, {θ, 0, 1}] // Normal // FullSimplify
Collect[% // Expand, Csc[_]]
(* (2^(-1 - m) θ^m C[
2] Csc[(l + m) π] Gamma[-m] Gamma[-l + m] Sin[(l - m) π])/Gamma[-l - m] + (
2^(-1 + m) θ^-m Gamma[m] (π C[2] Cos[m π] + 2 C[1] Sin[m π]))/π *)
(* 2^(-1 + m) θ^-m C[2] Cos[m π] Gamma[m] + (
2^(-1 - m) θ^m C[
2] Csc[(l + m) π] Gamma[-m] Gamma[-l + m] Sin[(l - m) π])/Gamma[-l - m] + (
2^m θ^-m C[1] Gamma[m] Sin[m π])/π *)
(* 总之我们由此分析出了C[2]必须为零。最后再反变换： *)
inverseFiniteFourierTransform[tsol /. C[2] -> 0, m, {φ, -Pi, Pi}] // Simplify

(* 这个结果相较于 SphericalHarmonicY 还差了系数，因为我们没归一化，不过LZ只是问这个方程怎么解，我就不继续弄了（别忘了特征函数乘以任意非零常系数后依旧是特征函数） *)
(*最后画个图：*)
Table[If[Abs@贴吧用户_00000eC🐾 <= l,
SphericalPlot3D[
E^(I m φ) LegendreP[l, m, Cos[θ]] // Abs, {θ, 0,
Pi}, {φ, 0, 2 Pi}],], {l, 0, 2}, {m, -2, 2}] // GraphicsGrid