I have the eigenvalues problem on $[0,\pi] \times [0,2\pi]$ $$\left(\frac{1}{\sin\theta}\frac{\partial}{\partial \theta} \left[\sin\theta \frac{\partial}{\partial \theta}\right] + \frac{1}{\sin^2\theta}\frac{\partial^2}{\partial \phi^2}\right)f (\theta,\phi) + g(\theta)f(\theta,\phi) = \lambda f(\theta,\phi).$$
The differential operator is the angular part of the laplacian in spherical coordinates and $g \in C^{\infty}$. What I want to know is: Are there any other techniques of solving this PDE than using the particular ansatz: $f(\theta,\phi) = h(\theta) \frac{e^{in \phi}}{\sqrt{\sin(\theta)}}$?
Also another ansatz would be interesting. I am just interested in any technique that could be useful here.
If anything is unclear, please let me know.
I would ordinarily leave a comment in this case, because I don't have any real answers. But it's too long to leave in a comment.
For ODEs $-ay''+by'+cy=\lambda w y$, the common techniques are
Multiply on the left by $\rho$ to get $\rho(ay''+by')= ((\rho a) y')'$ which (essentially) uniquely determines $\rho$. This requires minimal smoothness of coefficients.
Change the dependent variable $y=\rho f$ to obtain a new equation where $a$, $w$ are unchanged, but you end up with selfadjoint form $-(ay')'+qy=\lambda wy$. This one is not unique because you can multiply the entire equation by any factor before applying the technique, which is very helpful, but leaves you without unique decisions. This method also requires more smoothness for the coefficients, which limits its theoretical applications, but does not affect you for this problem.
Starting with $-(ay')'+qy=\lambda wy$ with $a > 0$, change the independent variable by writing $y(x)=f(\int \frac{1}{a}\,dx)$ assuming $1/a$ is locally integrable on the interior of the interval. Then you get $-\frac{1}{a}f''+qy=\lambda wy$, but you must invert $\int\frac{1}{a}\,dx$ to write this as $-f''+Qy=\lambda W y$.
Using these forms you can get $-f''+Qy=\lambda y$ with enough differentiability for the coefficients of the original equation.
About the only thing left to try is $f(\theta,\phi)=F(\int_{\pi/2}^{\theta}\frac{1}{\sin\theta'}\,d\theta',\phi)$ to transform your equation to $$ \frac{1}{\sin^{2}\theta}\left[F_{11}+F_{22}\right]+g(\theta)F=\lambda F. $$ You'll need to invert $u=\int_{\pi/2}^{\theta}\frac{1}{\sin\theta'}d\theta'$ to change $g(\theta)$ into $h(u)$, and $\sin^{2}(\theta)$ into $k(u)$. The new equation will be on $(-\infty,\infty)\times[0,2\pi]$. There may be something to be gained because of this becoming the Laplacian again, except that you'll have the new weighted space with weight $k(u)$ instead of $\sin\theta$, and some of the niceness of $g$ may be erased. I'm not hopeful because others have been down this road, and it sounds like you need something new.
Have you tried SLEIGN2 ? (Sturm-Liouville Eigenvalue/Eigenfunction package) This is a FORTRAN project by some impressive Mathematicians in the field. http://www.math.niu.edu/SL2/intro.pdf . Remove the intro.pdf from the link to get the *.f files.