\begin{equation} \hat{l}^2 Y(\theta , \varphi)=cY(\theta , \varphi) \end{equation}
Substituting the operator $\hat{l}^2$ into sphericals and writing the spherical harmonic as a product $S(\theta)T(\varphi)$ \begin{equation} -\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} \end{equation}
Acting with the operator on the function and simplifying:
\begin{equation} -\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}} \end{equation}
\begin{equation} \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) \end{equation}
Solving this differential equation gives us the eigenvalue, c, and the eigenfunction $S(\theta)$.
\begin{equation} c=\hbar^2l(l+1) \end{equation} \begin{equation} 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 \end{equation}