You can also do this with the numpy polynomial package.
In [1]: from numpy.polynomial import Polynomial, Legendre
In [2]: for i in range(5):
...: p = Legendre.basis(i).convert(kind=Polynomial)
...: print p.coef
...:
[ 1.]
[ 0. 1.]
[-0.5 0. 1.5]
[ 0. -1.5 0. 2.5]
[ 0.375 0. -3.75 0. 4.375]
Note that the coefficients go from low to high order. However, there is no need to do the conversion to a power series to compute the values and it is more accurate not to.
In [3]: Legendre.basis(2)(np.linspace(-1,1,10))
Out[3]:
array([ 1. , 0.40740741, -0.03703704, -0.33333333, -0.48148148,
-0.48148148, -0.33333333, -0.03703704, 0.40740741, 1. ])
You can also plot the result from [-1, 1] using the linspace method.
In [4]: plot(*Legendre.basis(2).linspace())
Out[4]: [<matplotlib.lines.Line2D at 0x30da4d0>]