Consider a Laplacian: $$\hat L=\frac1{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac1{\sin^2\theta}\frac{\partial^2}{\partial\varphi^2},$$ where $\theta\in[0,\pi], \varphi\in[0,2\pi)$.
I know that its eigenvalue problem (with boundary conditions as periodic in $\varphi$ and bounded in $\theta$) can be solved by separation of variables and the solution is spherical harmonics. Now however, I'm interested in treating this problem numerically (because I need to later be able to solve eigenvalue problem for $\hat L+U(\theta,\varphi)$).
I want to use a grid method like finite differences for this.
What I found is that for $\theta\to0$ and $\theta\to\pi$ coefficients go to infinity as $\cot\theta$ for first part and $\csc^2\theta$ for second one. So, using a grid method, I have to reduce $\theta$ domain to some $[\varepsilon,\pi-\varepsilon]$. This of course still leads to bad representation of $\cot$ and $\csc^2$ at ends. So, I try to substitute $\cot\theta=z$, then $\csc^2\theta=z^2+1$ to have more fine discretization near poles.
Now new Laplacian looks like $$\hat {L^\prime}=(z^2+1)^2\frac{\partial^2}{\partial z^2}+z(z^2+1)\frac{\partial}{\partial z}+(z^2+1)\frac{\partial^2}{\partial\varphi^2}.$$ Here $z\in(-\infty,\infty)$; for numerical treatment I replace infinities with some $z_{max}$. Now I mostly successfully solve this problem, BUT: it appears that although it's enough for $\varphi$ to have 10 points on grid, for $z$ even 300 points gives too low precision, whatever $z_{max}$ I select (though of course, getting them higher gives better accuracy, and $z_{max}$ has to be selected not too big and not too small to discretize space near poles and far from them equally well).
At the same time, the exact eigenfunctions look really simple and looking at them I wouldn't even think there's some difficulty at poles. This makes me think there must be much simpler way to solve this problem than using that many points in $\theta$.
So, my question: is there a good way to transform $\hat L$ so that discretization with small number of points in $\theta$ would give good accuracy - similar to what I have for $\varphi$ now?
