Solve the angular part of Schrodinger equation numerically I would like to solve the angular part (the one for what is usually called the $\theta$ angle) of a time-independent 3D Schrodinger equation 
$$
\frac{\mathrm{d}}{\mathrm{d}x}\left[ (1-x^2) \frac{\mathrm{d}P(x)}{\mathrm{d}x}
\right]+\left[ l(l+1) - \frac{m^2}{1-x^2} \right]P(x) = 0,
$$
where $l=0,1,2,\ldots$ and $m = -l, -l+1, \ldots, l$ as usual and $x\in[-1,1]$,
Now, the complication is that I want to do it numerically. Analytically, one
gets a bunch of Legendre polynomials and spherical harmonics. However, for me it is unclear which boundary conditions should I set. 
One boundary condition will probably be equivalent to the normalization of my solutions. In order to make it compatible with the Legendre polynomials, I can set it to
$$ P(1) =1. $$
However, what about the second one (it is a second-order ODE after all)? I guess, it should somehow encode the fact that my solutions should be bounded.
Any comments, including sending me to RTFM (with appropriate links) are more than welcome!
 A: What you are doing is an eigenvalue problem. Eigenvectors are determined by the space you are looking at, and this is why you usually specify some boundary conditions. In your case just the requirement of absence of singularities should do the job (i.e. you want some subspace of $C[-1,1]$). This is the analytical viewpoint. 
The numerical viewpoint actually depends on your algorithm. First of all, if you really want "to solve the equation numerically", I assume that you are playing the game of not knowing the answer. So you do not actually know that $l$ is integer beforehand. If I were solving the problem, I would put it on a lattice and then write it as a finite-dimensional eigenvalue problem. In deriving the finite difference equations I would use the fact that my solution should be finite at the endpoints of the interval.
A way to do this is to introduce homogenious lattice at points (lets call them so) 0,1,2,3,4... Then integrate the equation from $i-1/2$ to $i+1/2$ and use middle-difference fromulae for derivatives (you will need their values at $i\pm1/2$, so the middle difference will return you back to your lattice) and middle-rectangles formula for integrals (it is important to use approximations of the proper order. I believe that I am telling you an algorithm of second-order presicion). Then you will have to do something with the endpoints. For them do the integration from $0$ to $1/2$ and respectively on the other end. In doing so you will use the fact that $(1-x^2)\frac{dP}{dx}$ is $0$ at the endpoints. And this picks up the appropriate space for your solutions.
Long story made short, I believe that at least for some calculational schemes the conditions should be that $(1-x^2)\frac{dP}{dx}$ vanishes at the endpoints.
A: Since the coordinate is an angle, you should specify the periodic boundary conditions.
