I have a physical problem with a solution component that involves an infinite series of Legendre polynomials. I'm trying to compute the coefficients but am getting some odd behaviour. As I understand it, provided $-1 \leq x \leq x$, then for $f(x)$ I can get my coefficients by
$$A_{n} = \frac{2n+1}{2} \int_{-1}^{1} P_{n}(x)f(x)\space dx.$$
where $P_{n}$ is the legendre function of degree $n$. My functions $f(x)$ uses cosine functions (which satisfy $|x| \leq 1$). Writing $x = \cos\theta$, then an example of the type of function we need to solve is
$$f(x) = \frac{100}{\sqrt{25 + x^2}}.$$
However, when I compute these, I get some odd behavior - for all odd values, $A_{n} = 0$, which is fine. Intially, it seems like the series converges (as expected) with the first few non-zero terms gradually getting smaller and alternating sign;
$A_{2} = -0.26$, $A_{4} = 0.26$, $A_{6} = -2.596 \times 10^{-4}$, $A_{8} = -1.266 \times 10^{-6}$....
But this does not continue -after this, the values rapidly increase, as shown in the figure below, eventually exploding to huge values as $n \rightarrow \infty$ 
My question is why does this occur, and where am I going wrong? I've seen conflicting things about whether $f((x)$ has to be a polynomial for this to work, so if that's the problem, how can I express the function $f(x)$ if it is non-polynomial, or is there something more subtle (or obvious) I'm missing?
So part of the problem was apparently numerical instability - changing the language I used and the errors bounds made the function behave much more predictably. It seems the integrals were blowing up in the original code, which disappeared when they were solved in Mathematica. Though if there are any conditions I should know about for $f(x)$, I'd still be grateful to hear them!