$$\hat{l}^2 Y(\theta , \varphi)=cY(\theta , \varphi)$$

Substituting the operator $\hat{l}^2$ into sphericals and writing the spherical harmonic as a product $S(\theta)T(\varphi)$ $$-\hbar^2\left(\frac{\partial^2}{\partial\theta^2}+cotg\theta \frac{\partial}{\partial\theta}+\frac{1}{sin^2\theta}\frac{\partial^2}{\partial\varphi^2}\right)\left(S(\ theta)\frac{1}{\sqrt{2\pi}}e^{im\varphi}\right)=cS(\theta)\frac{1}{\sqrt{2\pi}}e^{im \varphi}$$

Acting with the operator on the function and simplifying:

$$-\hbar^2\left(\frac{\partial^2S(\theta)}{\partial\theta^2 }+cotg\theta\frac{\partial S}{\partial\theta}-\frac{m^2}{sin^2\theta}S(\theta)\right)\cancel{\frac{1}{ \sqrt{2\pi}}e^{im\varphi}}=cS(\theta)\cancel{\frac{1}{\sqrt{2\pi}}e^{im\varphi}}$$

$$\frac{\partial^2S(\theta)}{\partial\theta^2}+cotg\theta\frac{\partial S}{\partial\theta}-\frac{m^ 2}{sin^2\theta}S(\theta)=-\frac{c}{\hbar^2}S(\theta)$$

Solving this differential equation gives us the eigenvalue, c, and the eigenfunction $S(\theta)$.

$$c=\hbar^2l(l+1)$$ $$S_{l,m}(\theta)=(-1)^{\frac{m+|m| }{2}}sin^{|m|}\theta\sum_{j=0,(1),2,(3)..}^{l-|m|}a_j cos^{j}\theta$$