if you have a free Schrödinger operator on a sphere $-\Delta \psi(\theta,\phi) = E\psi(\theta,\phi),$ then the substitution $\psi(\theta,\phi) = f(\theta)e^{i n \phi}$ leaves you with the differential equation
$$- f''(\theta)- \cot(\theta) f'(\theta) + \frac{n^2}{\sin^2(\theta)} f(\theta)+ = E f(\theta).$$
For $n=0$ this simplifies to $$- f''(\theta)- \cot(\theta) f'(\theta)= E f(\theta).$$
My question is: What are the correct boundary conditions for the ODE after the substitution?
The Calculus distance scale factors in spherical coordinates are $s_{r}=1$, $s_{\theta}=r$, $s_{\phi}=r\sin\theta$. So the Laplacian is $$ \Delta=\frac{1}{s_{r}s_{\theta}s_{\phi}}\left[ \frac{\partial}{\partial r}\frac{s_{\theta}s_{\phi}}{s_{r}}\frac{\partial}{\partial r} +\frac{\partial}{\partial\theta}\frac{s_{r}s_{\phi}}{s_{\theta}}\frac{\partial}{\partial\theta}+\frac{\partial}{\partial\phi}\frac{s_{r}s_{\theta}}{s_{\phi}}\frac{\partial}{\partial\phi}\right]\\ %% = \frac{1}{r^{2}\sin\theta}\left[\frac{\partial}{\partial r}r^{2}\sin\theta\frac{\partial}{\partial r} %% + \frac{\partial}{\partial\theta}\sin\theta\frac{\partial}{\partial\theta} %% + \frac{\partial}{\partial\phi}\frac{1}{\sin\theta}\frac{\partial}{\partial\phi}\right]\\ = \frac{1}{r^{2}}\frac{\partial}{\partial r}r^{2}\frac{\partial}{\partial r}+ \frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\sin\theta\frac{\partial}{\partial\theta}+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}. $$ The following restriction to the sphere is the correct Laplacian on the sphere: $$ \Delta_{S}=\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\sin\theta\frac{\partial}{\partial\theta}+\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}. $$ The natural $L^{2}_{S}$ inner-product is $$ (f,g)_{S} = \frac{1}{4\pi}\int_{0}^{\pi}\int_{0}^{2\pi}f(\theta,\phi)\overline{g(\theta,\phi)}\sin\theta\,d\phi d\theta. $$ And, as you noted, if $\psi(\theta,\phi)=\Theta(\theta)e^{in\phi}$, then $$ \{-\Delta+V(\theta)\}\psi = e^{in\phi}\left[-\frac{1}{\sin\theta}\frac{d}{d\theta}\sin\theta\frac{d}{d\theta}-\frac{n^{2}}{\sin^{2}\theta}+V(\theta)\right]\Theta(\theta). $$ The function $\phi(\theta,\phi)=e^{in\phi}\Theta(\theta)$ is in $L^{2}_{S}$ iff $\Theta \in L^{2}_{\sin}[0,\pi]$, which is to say $$ \|\Theta\|_{\sin}^{2}=\int_{0}^{\pi}|\Theta(\theta)|^{2}\sin\theta\,d\theta < \infty. $$ It is not enough to require that $\psi \in L^{2}_{S}$ in order to study this problem in a classical manner. For example, suppose $V=0$ and suppose $n=0$; then $$ \Delta \Theta(\theta) = \frac{1}{\sin\theta}\frac{d}{d\theta}\sin\theta\frac{d}{d\theta}\Theta. $$ Therefore $\Delta\theta(\theta)=0$ iff $$ \Theta(\theta) = A+B\ln\left(\frac{1-\cos\theta}{1+\cos\theta}\right). $$ Note: All such solutions are in $L^{2}_{S}$, but you don't want the solution where $A=0$, $B=1$ being part of a classical solution. The presence of the logarithmic solution is an artifact of the vanishing Jacobian $\sin\theta$ of the coordinate transformation; this vanishing allows the fundamental solution of $\Delta F=\delta_{n}-\delta_{s}$ to creep into things, where $\delta_{n}$ is the Dirac delta at the North pole and $\delta_{s}$ is the Dirac delta at the South pole.
So why is it not a good idea to allow the logarithmic solution into your formulation of the problem? Physicists sometimes argue that it makes no sense to allow an unbounded pure state in the system. There's a compelling case to made for this, but it's not entirely satisfying to invoke such a condition when there is a clear Mathematical necessity to disallow this type of state based on the axioms of Quantum Mechanics: If you define $\Delta$ on the domain consisting of $C^{2}$ functions on the manifold, and you add the logarithmic function to the domain, then the resulting operator is no longer symmetric! And that leads to an unacceptable Hamiltonian for Physics. Indeed, there are $C^{2}$ functions $\psi$ on $\mathbb{S}^{2}$ such that $$ \left(\Delta\ln\left(\frac{1+\cos\theta}{1-\cos\theta}\right),\psi\right)_{S}=0 \ne \left(\ln\left(\frac{1+\cos\theta}{1-\cos\theta}\right),\Delta\psi\right)_{S}. $$ That is reason enough to exclude the logarithmic function from any domain of a Hamilton which includes $C^{2}$ functions on the manifold: Quantum Hamiltonians are selfadjoint and, hence, also symmetric. The exclusion of this function is equivalent to the requirement that pure states be bounded, which justifies the Physicist's reasoning.
A standard way to arrive at a selfadjoint operator $\Delta_{S}$ is to start with functions that are $C^{2}$ on the sphere as a manifold, and then show that $\Delta$ is symmetric and has a sefladjoint closure $\Delta_{S}$; that's one way to overcome the artifact of a vanishing Jacobian. And it leads to an alternative characterization of the domain of $\Delta_{S}$ as consisting of functions in a Sobelev space. This is because the first order derivatives remain in $L^{2}_{S}$ under this closure operation because of the identity $$ \begin{align} (-\Delta f,f)_{S} & = \int_{S}|\nabla f|^{2}\,dS \\ & = \int_{S}\left( \left|\frac{\partial f}{\partial\theta}\right|^{2} +\left|\frac{1}{\sin\theta}\frac{\partial f}{\partial\phi}\right|^{2}\right)\,dS < \infty. \end{align} $$ And this condition excludes the fundamental logarithmic solution because $$ \frac{\partial}{\partial\theta}\ln\left(\frac{1-\cos\theta}{1+\cos\theta}\right)=\frac{1}{\sin\theta} \notin L^{2}_{S}. $$ The resulting operator is unique inasmuch as it is the closure from a natural manifold domain which is invariant under all symmetries of $\mathbb{S}^{2}$. This one gradient condition is then propogated to the individual eigenfunction equations: $$ \int_{0}^{\pi}\left(|\Theta_{n,m}'(\theta)|^{2}+\frac{n^{2}}{\sin^{2}\theta}|\Theta_{n,m}(\theta)|^{2}\right)\sin\theta\,d\theta < \infty. $$ This last condition is met for all of the classical Associated Legendre Polynomial solutions $\Theta_{n,m}$, in addition to the requirements that $\Theta_{n,m}$ be in $L^{2}_{\sin}[0,\pi]$.
Added: Back to your case. If $V$ is smooth in $\theta$, then you can view the addition of a potential as a bounded perturbation with $\|V\|\le \sup|V|$. This is a general fact: If $A$ is a bounded on unbounded selfadjoint operator and $V$ is a bounded selfadjoint operator then $A+V$ is selfadjoint on $\mathcal{D}(A)$. So I'm not going to worry about $V$ in the discussion. Without knowing more about your specific $V$--which has not been answered--I was prepared to deal with forms in order to handle the more general perturbation.
As a final note, I'll be more explicit about the separation of variables connection. You can write every $f \in L^{2}_{S}$ as $\sum_{n=-\infty}^{\infty}f_{n}(\theta)e^{in\phi}$, and $\|f\|^{2}_{S}=\sum_{n=-\infty}^{\infty}\|f_{n}\|_{L^{2}_{\sin}[0,\pi]}^{2}$ follows.
To put this back on the ODEs, let $-\Delta_{0}$ be the restriction of $-\Delta$ to the $C^{\infty}$ functions which vanish in a neighborhood of the North and South poles. If $f$ is $\mathcal{D}(-\Delta_{0})$, then $f=\sum_{n=-\infty}^{\infty}e^{in\phi}f_{n}(\theta)$, where $f_{n}$ is $C_{c}^{\infty}(0,\pi)$. For any such $f \in \mathcal{D}(-\Delta_{0})$, $$ -\Delta_{0} f = \sum_{n=-\infty}^{\infty}e^{in\phi}L_{n}f_{n},\\ L_{n}f_{n}=\left[-\frac{1}{\sin\theta}\frac{d}{d\theta}\sin\theta\frac{d}{d\theta}+\frac{n^{2}}{\sin^{2}\theta}\right]f_{n}, $$ and $$ \begin{align} (-\Delta_{0} f,f)_{S} & = \sum_{n=-\infty}^{\infty} (L_{n}f_{n},f_{n})_{\sin} \\ & = \sum_{n=-\infty}^{\infty} \|f_{n}'\|_{\sin}^{2}+\|\frac{n}{\sin\theta}f_{n}\|^{2}_{\phi}. \end{align} $$ The Friedrichs extension $-\Delta_{f}$ of $-\Delta_{0}$ is the restriction of $(-\Delta_{0})^{\star}$ to the domain of functions $f$ for which the above form is finite. The Friedrichs extension $-\Delta_{f}$ is automatically selfadjoint and semibounded.
The nice thing about the Friedrichs extension is that the conditions for it can be put off onto the ordinary differential operators (ODOs) through the individual quadratic forms, which is evident from the sum representation for $(-\Delta_{0}f,f)_{S}$. This gives us something substantial because the adjoint of an ODO is well known and classical, whereas the domain for the PDO is not so easy to study in general. For example, the Friedrichs extension of $L_{n}$ has a domain consisting of all twice locally absolutely continuous functions $f$ for which $f, L_{n}f, f', \frac{n}{\sin}f$ are in $L^{2}_{\sin}$. On this domain integration by parts is allowed, and the conditions can be translated to endpoint conditions. For example, suppose $n\ne 0$. Then $$ f_{n}, \frac{1}{\sin\theta}f_{n}', -\frac{1}{\sin\theta}(\sin\theta f_{n}')'+\frac{n}{\sin^{2}\theta}f_{n} \in L^{2}_{\sin}. $$ So all of these, and the products of any two of these are in $L^{1}_{\sin}[0,\pi]$. That gives you a lot of things to look at, and they all result in a limiting expressing for an integral with respect to the upper limit. For example, the following has finite limits as $\theta$ approaches either endpoint: $$ f_{n}(\theta)-f_{n}(\pi/2)=\int_{\pi/2}^{\theta}\frac{1}{\sin\theta'}f_{n}'(\theta')\,\sin\theta'\,d\theta' $$ A condition more directly related to the form is $$ \int_{\pi/2}^{\theta} \overline{f_{n}}L_{n}f_{n}\,d\theta' = \int_{\pi/2}^{\theta}\left( -\overline{f_{n}(\theta')}(\sin\theta'f_{n}'(\theta'))'+\frac{n}{\sin\theta'}|f_{n}|^{2}\right)d\theta' \\ = \left.-\sin\theta' \overline{f_{n}}f_{n}'\right|_{\pi/2}^{\theta} + \int_{\pi/2}^{\theta}\sin\theta'|f_{n}'|^{2}+\frac{n}{\sin\theta}|f_{n}|^{2}\,d\theta' $$ From this last condition you can conclude that the evaluation term has a limit as $\theta$ approaches either endpoint. Similarly, the following also has endpoint limits: $$ \sin\theta \frac{d}{d\theta}|f_{n}|^{2}. $$ And, the endpoint limits are $0$. Otherwise there is a non-zero $C$ such that $$ \frac{d}{d\theta}|f_{n}|^{2} \approx \frac{C}{\sin\theta}, $$ which will contradict the fact that $|f_{n}|^{2}$ has a finite limit at each endpoint. I'll let you continue the analysis to match Teschl's.
Discrete spectrum of $L_{n}+V$: The spectrum of $L_{n}$ is discrete with eigenvalues $\lambda_{n,m}$ that tend to $\infty$. If $\lambda\notin\sigma(L_{n})$, then $(L_{n}-\lambda I)^{-1}$ is compact because it has a complete orthonormal basis of eigenfunctions with eigenvalues $\frac{1}{\lambda_{n,m}-\lambda}$, a sequence which tends to $0$ as $m\rightarrow\infty$. The resolvent of $L_{n}+V$ is equal to the following for all $\lambda$ for which $\|V(L_{n}-\lambda I)^{-1}\|< 1$: $$ \begin{align} (L_{n}+V-\lambda I)^{-1} & =\left[(L_{n}-\lambda I)(I+(L_{n}-\lambda I)^{-1}V)\right]^{-1} \\ & = (I+(L_{n}-\lambda I)V)^{-1}(L_{n}-\lambda I)^{-1} \end{align} $$ In particular, because $\|(L_{n}-\lambda I)^{-1}\|\le 1/|\Im\lambda|$, then the above holds at least for $\|V\|_{\infty} < |\Im\lambda|$. That's enough to conclude that some resolvent operator $(L_{n}+V-\lambda I)^{-1}$ is compact, which guarantees the discreteness of the spectrum of $L_{n}+V$.