Hi, I have a symbolic partial differential equation like $\nabla^2 P(t,r,\theta,\phi) = 0$. Now I want to substitute the following anzast $P(t,r,\theta,\phi) = e^{-i\omega t} f(r) Y^{m}_{l}(\theta,\phi)$ where the last terms is spherical bessel function. When I wrote,

```
%display latex
M = Manifold(4, 'M', structure='Lorentzian')
X.<t,r,th,phi> = M.chart(r"t r:(0,+oo) th:(0,pi):\theta phi:(0,2*pi):\phi ")
g = M.metric()
g[0,0] = 1-2/r
g[1,1] = -1/(1-2/r)
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
w = var('w')
T = e^(-i*w*t)
f = M.scalar_field(function('f')(r), name='f')
l = var('l')
m = var('m')
Y = spherical_harmonic(l,m,th,phi)
P = T*f*Y
eq1 = (dalembertian(P).expr()) == 0
eq1.simplify()
```

I should get an expression independent of $\theta,\phi$ since angular part of the equation cancel itself because by definition spherical harmonics satisfies it's own differential equation (see Wolfram-Spherical Harmonics for example). Now when I check manually using

```
Q = (1/(sin(th))*diff(sin(th)*diff(Y,th),th)+1/((sin(th))^2)*diff(diff(Y,phi),phi)+l*(l+1)*Y)
bool(Q==0)
```

I got FALSE. To force Sage to use this relation, I tried things like,

```
eq1.substitute_expression(Q==0).simplify_full()
```

but they did not work either. The problem roots from the following line,

```
differentiate(spherical_harmonic(l,m,th,phi),th)
```

gives an expression in terms of higher order harmonics and Sage is unable to simplify expressions using recursion relations of the spherical harmonics. I also looked at the following integrals

```
integrate(integrate(spherical_harmonic(l,m,th,phi)*conjugate(spherical_harmonic(l,m,th,phi))*sin(th),th,0,pi),phi,0,2*pi).simplify_full()
```

but again Sage is unable to evaluate the types of integral unless we plug specific $l,m$. My question is there any way to defince spherical harmonics without using built-in spherical_harmonic() function (like power series) so that Sage would be able to simplify "it's differential equation" involving spherical harmonics to automatically zero.