I have studied Theoretical Chemistry and Physics. While I was doing research on numerical integration I found out that zeroes of orthogonal polynomials are important. Mathematica became general available around that time at universities and I decided to take a look at the distribution of the zeroes of Legendre polynomials so I typed:

This plots the N zeroes of the N-th Legendre polynomial as points (x, y), where x is the index 1..N scaled back to the interval 0..1, and y is the zero for the first 20 Legendre polynomials. As you can imagine I was quite surprised to see the following plot:

It clearly illustrates the power of mathematical visualisation. Since the zeroes are distributed as a sine function, the ArcSin of the zeroes should almost lie on a straight line and sure enough:

So a first order approximation to the zeroes is Sin(a*x+b). Since the zeroes of Legendre polynomials are symmetric around 0, the line should be symmetric around x=0.5, so the first order approximation can be simplified to Sin(a*(x-1/2)). Let’s see how the coefficient a that determines the line depends on N. Redoing the calculation for the Legendre polynomials 2..60, and now for each polynomial fitting the ArcSin of the zeroes to a line of the form a*(x-1/2) gives:

The coefficient a seems to depend in a well-behaved way on N, asymptotically approaching the value Pi. It looks like an inverted and shifted ArcTan, so the Tan of -a-Pi/2 might look like a straight line, and sure enough:

A better approximation for the ArcSin of the i-th zero of the N-th Legendre polynomial finally becomes:

(-ArcTan(0.5748+N*0.6442) + 3*Pi/2)*(i/(N+1)-1/2)

This approximation is accurate to within 0.01 for the first 60 Legendre polynomials. In fact, the largest deviations occur for the lowest order polynomials 2, 3 and 4. From order 5 and up the approximation is accurate to within 0.002 as the following figure shows: