I'm writing a program that solves an integral by Legendre-Gauss quadrature. The algorithm for nth-order quadrature requires, at one point, finding the roots of the nth-order Legendre polynomial, Pn(x), assigning them to the array Absc (for 'abscissa'). Pn is an nth order polynomial with n independent real roots on the interval [-1,1]. I'd like to be able to compute the roots, instead of just importing them from some library. I am able to create an array that gives the polynomial coefficients, which I call PCoeff. To find the roots I have tried

```
 Absc = numpy.roots(PCoeff)

```

This works up to about n = 40, but beyond that it starts to fail, giving complex roots when it really shouldn't be. I've also tried defining the polynomial using

```
P = numpy.poly1d(PCoeff)
Absc = P.r

```

but this gives the same problems, presumably because it uses the same numpy root-finding algorithm.

Another method, which seemed promising was using scipy.optimize.fsolve(Pn, x0), where x0 is an n-element array of my guess at the roots. The problem with this is that, depending on my x0 choices, this method may give one particular root more than once in the place of other roots. I've tried populating x0 as equidistant points on [-1,1]

```
x0 = numpy.zeros(n)
step = 2./(n-1)
for i in xrange(n):
  x0[i] = -1. + i*step

```

but once I get to n = 5, fsolve gives some roots repeated and neglects others. I have also tried using the results of numpy.roots as x0. However, in the problem area where np.roots gives complex values, these cause an error in fsolve

```
TypeError: array cannot be safely cast to required type

```

I saw online that there is a scipy.optimize.roots() routine that could work, but it is not in the scipy library on my computer. Updating is a hassle since I don't have permission to download things on this computer.

I'd like to be able to run the quadrature at 64th order for high accuracy, but this root finding causes failure. Any ideas?