Given the recursively defined sequence $$ a_2 = 2(C+1)a_0 $$ and $$ a_{n+2}=\frac{\left(n(n+6)+4(C+1)\right)a_n - 8(n-2)a_{n-2}}{(n+1)(n+2)}. $$ that I got from the Frobenius method applied to an ODE, where for every odd $n$ we have $a_n=0$. This brought me to the following problem:
For which $C\in \mathbb{R}$ do we have that $|\sum_{n=0}^{\infty}a_n |< \infty$? Of course, there is a trivial solution $C=-1$, but are there more available?
Edit: Furthermore, it was explained in the comments below that the convergence (due to linearity) is independent of the initial value $a_0$, so we could also just set this one $=1$.
If anything is unclear, please let me know. Since nobody answered so far, I just wanted to say that also interesting observations are highly appreciated.
Your reduced equation is $$ ( -x^2+1 ) \frac {{\rm d}^2}{{\rm d}x^2} f ( x ) -x \frac {\rm d}{{\rm d}x} f ( x ) +4 ( 2x^2-1+ ( 2x^2-1 )^2+C ) f ( x ) =0. $$ Symmetry and Spectral Theory: This is not in symmetric Sturm-Liouville form, which is required in order to do meaningful spectral theory. There are two ways to put it into symmetric form: by letting $f=\rho g$ and choosing $\rho$, or by multiplying everything on the left by some $\rho$. The second is easier, but the first method is in keeping with how to transform the associated Legendre equation into the form encountered through separation of variables of the Laplacian in spherical coordinates. So, that's the one I'll use, because of this parallel. Let $f=\rho g$, and divide by $\rho$ on the left. The above equation becomes $$ \begin{align} &(1-x^{2})g''+2(1-x^{2})\frac{\rho'}{\rho}g'+(1-x^{2})\frac{\rho''}{\rho} g \\ &-xg'-x\frac{\rho'}{\rho}g + \\ &+4(2x^{2}-1+(2x^{2}-1)^{2}+C)g=0. \end{align} $$ Now choose $\rho$ so that the terms involving $g''$, $g'$ become $\frac{d}{dx}(1-x^{2})\frac{d}{dx}$. Equivalently, $$ 2(1-x^{2})\frac{\rho'}{\rho}-x=-2x,\\ \frac{\rho'}{\rho}=\frac{-x}{2(1-x^{2})}=\frac{1}{4}\frac{(1-x^{2})'}{(1-x^{2})},\\ \rho = (1-x^{2})^{1/4}. $$
Then the equation becomes $$ \frac{d}{dx}(1-x^{2})\frac{d}{dx}g+(1-x^{2})\frac{\rho''}{\rho}g-x\frac{\rho'}{\rho}g+4(2x^{2}-1+(2x^{2}-1)^{2}+C)g=0 $$ We already have $\rho'/\rho$. Differentiating $\rho'/\rho$ gives $$ \begin{align} \left(\frac{\rho'}{\rho}\right)' & =\frac{\rho''}{\rho}-\left(\frac{\rho'}{\rho}\right)^{2}, \\ \frac{\rho''}{\rho} & =\left(\frac{\rho'}{\rho}\right)^{2}+\left(\frac{\rho'}{\rho}\right)' \\ & = \frac{x^{2}}{4(1-x^{2})^{2}}-\frac{2(1+x^{2})}{4(1-x^{2})^{2}} = -\frac{2+x^{2}}{4(1-x^{2})^{2}} \end{align} $$ Finally, $$ (1-x^{2})\frac{\rho''}{\rho}-x\frac{\rho'}{\rho}=-\frac{2+x^{2}}{4(1-x^{2})}+\frac{x^{2}}{2(1-x^{2})}=-\frac{2-x^{2}}{1-x^{2}}=-1-\frac{1}{1-x^{2}}. $$ The transformed equation becomes $$ \frac{d}{dx}(1-x^{2})\frac{d}{dx}g-\left(1+\frac{1}{1-x^{2}}\right)g+4(2x^{2}-1+(2x^{2}-1)^{2}+C)g=0. $$ The Legendre equation obtained from the Laplacian through separation of variables in spherical coordinates looks like $$ \left[\frac{d}{dx}(1-x^{2})\frac{d}{dx}-\frac{m^{2}}{1-x^{2}}+n(n+1)\right]g = 0. $$ So this resembles the Associated Legendre Equation of order $m=1$. And the solutions you want will have the form $f=\rho g=(1-x^{2})^{1/4}g$. It is reasonable to expect that the values of $C$ should be asymptotic to $\frac{1}{4}(n(n+1)-1)$. Solutions $g_{1}$, $g_{2}$ for distinct values of $C$ will be orthogonal in the sense that $$ \int_{-1}^{1}g_{1}g_{2}\,dx = 0, $$ provided that you can find $g_{j}$ which are square-integrable at both endpoints (I suspect you can because of the symmetry you noted.) So, in the setting in which you are interested, it is natural to seek solutions $f_{j}=(1-x^{2})^{1/4}g_{j}$ which satisfy $$ \int_{-1}^{1}(1-x^{2})^{-1/2}f_{j}^{2}\,dx < \infty, $$ and you should expect to automatically have orthogonality in the following sense: $$ \int_{-1}^{1}(1-x^{2})^{-1/2}f_{1}f_{2}\,dx = 0. $$ Look into the Associated Legendre Polynomials http://en.wikipedia.org/wiki/Associated_Legendre_polynomials for $m=1$: $$ P_{n}^{m}=\frac{(-1)^{m}}{2^{n}n!}(1-x^{2})^{m/2}\frac{d^{m+n}}{dx^{m+n}}(x^{2}-1)^{n}. $$ These are the asymptotics I would expect: $f = (1-x^{2})^{1/4}g \sim (1-x^{2})^{1/4}P_{n}^{1}$. Note that for $m=1$, $P_{n}^{m}$ are not really polynomials.
Power Series Solutions: The goal of putting the equation into symmetric form, while it does lead to orthogonality and spectral (eigenvalue) information, seems to work against power series solutions in this particular case. To get rid of the singular terms at $x=\pm 1$, look instead at the indicial equation at these endpoints. It is enough to look at $x=-1$ because of symmetry. The equation has the approximate form $$ (1+x)(2)\frac{d^{2}}{dx^{2}}f+\frac{d}{dx}f=0. $$ The indicial equation is found by setting $f=(1+x)^{\rho}$ and solving for $\rho$ such that $$ 2\rho(\rho-1)+\rho=0\implies \rho=0 \mbox{ or } \rho=1/2. $$ The one with the larger real part is significant. So we expect that $f=(1-x^{2})^{1/2}g$ is the only form that could possibly lead to a power series $g$ which would be holomorphic across $x=\pm 1$. Here the goals of spectral theory and of finding nice solutions near the endpoints do not coincide. But, the same is true for the general Legendre equation, too. The functions $P_{n}^{1}$ have a factor of $(1-x^{2})^{1/2}$ times a polynomial. That's a sanity check that setting $f=(1-x^{2})^{1/2}g$ has a chance of leading to a polynomial $g$, an entire function $g$, or at least to function $g$ which remains controlled near the endpoints. Any such case would be significant to what you want, and the special values of $C$ required to gain control over $g$ are the expected eigenvalues.
So the next step is to substitute $f=(1-x^{2})^{1/2}g$ into the equation at the very top of this post. What I found is $$ \frac{1}{\sqrt{1-x^{2}}}\left((1-x^{2})\frac{d^{2}}{dx^{2}}-x\frac{d}{dx}\right)(\sqrt{1-x^{2}}g)=\left((1-x^{2})\frac{d^{2}}{dx^{2}}-3x\frac{d}{dx}-1\right)g. $$ The cancellation of negative powers of $1-x^{2}$ is nice and somewhat expected. The new equation for $g=f/\sqrt{1-x^{2}}$ is $$ (1-x^{2})g''-3xg'+4\left\{(2x^{2}-1)+(2x^{2}-1)^{2}+C-\frac{1}{4}\right\}g = 0,\\ (1-x^{2})g''-3xg'+(16x^{4}-8x^{2}+4C-1)g=0. $$ Set $g=\sum_{n=0}^{\infty}a_{n}x^{n}$ and substitute into the equation: $$ \sum_{n=0}^{\infty}(n+1)(n+2)a_{n+2}x^{n}-\sum_{n=0}^{\infty}n(n-1)a_{n}x^{n} -\sum_{n=0}^{\infty}3na_{n}x^{n}\\ +\sum_{n=4}^{\infty}16a_{n-4}x^{n}-\sum_{n=2}^{\infty}8a_{n-2}x^{n}+(4C-1)\sum_{n=0}^{\infty}a_{n}x^{n}=0. $$ The complete set of coefficient conditions leaves $a_{0}$ and $a_{1}$ arbitrary: $$ \begin{align} x^{0}:\; & 2a_{2}+(4C-1)a_{0}=0 \\ & a_{2}=\frac{1}{2}(1-4C)a_{0} \\ x^{1}:\; & 6a_{3}-3a_{1}+(4C-1)a_{1}=0 \\ & a_{3}=\frac{1}{6}(4-4C)a_{1} \\ x^{2}:\; & 12a_{4}-2a_{2}-6a_{2}-8a_{0}+(4C-1)a_{2}=0\\ & a_{4}=\frac{1}{12}(9-4C)a_{2}+\frac{8}{12}a_{0} = \frac{1}{4!}(1^{2}-4C)(3^{2}-4C)a_{0}+\frac{8}{(4)(3)}a_{0}\\ x^{3}:\; & 20a_{5}-6a_{3}-9a_{3}-8a_{1}+(4C-1)a_{3}=0\\ & a_{5}=\frac{1}{20}(16-4C)a_{3}+\frac{8}{20}a_{1} = \frac{1}{5!}(2^{2}-4C)(4^{2}-4C)a_{1}+\frac{8}{(5)(4)}a_{1}\\ x^{n}:\; & (n+1)(n+2)a_{n+2}-n(n-1)a_{n}-3na_{n}+16a_{n-4}-8a_{n-2}+(4C-1)a_{n}=0\\ & a_{n+2}=\frac{(n+1)^{2}-4C}{(n+1)(n+2)}a_{n}-\frac{8}{(n+1)(n+2)}a_{n-2} +\frac{16}{(n+1)(n+2)}a_{n-4}. \end{align} $$ I don't have enough time to push this forward to a pattern for the general term, but it looks promising ... will come back to it. This looks more like the Legendre arguments that I thought it would. Superficially it looks like the convergence is changed when $C$ takes on certain values. Please check if everything looks okay to you ...
Eigenfunction Conditions: This is a little messy. The original $f$ in your reduced equation gives rise to $g=(1-x^{2})^{-1/4}f$ which is a solution of $$ Lg=-\frac{d}{dx}(1-x^{2})\frac{d}{dx}g+\frac{1}{1-x^{2}}g-(16x^{2}-8x^{2}-1)g=4Cg. $$ The natural domain is where $(Lg,g)$ exists, where this is the $L^{2}[-1,1]$ inner-product. Integrating by parts gives $$ (Lg,g)=\|\sqrt{1-x^{2}}g'\|^{2}+\|g/\sqrt{1-x^{2}}\|^{2}+\cdots. $$ In other words, the natural conditions on the original $f$ are that $$ (1-x^{2})^{1/2}g'=(1-x^{2})^{1/2}((1-x^{2})^{-1/4}f)' \in L^{2}[-1,1],\\ (1-x^{2})^{-1/2}g=(1-x^{2})^{-3/4}f \in L^{2}[-1,1]. $$ The power series $p$ for which the coefficient equations are given above is $p=(1-x^{2})^{-1/2}f$. Or, $f=(1-x^{2})^{1/2}p$, so the $p$ should satisfy $$ (1-x^{2})^{1/2}((1-x^{2})^{1/4}p)'=(1-x^{2})^{3/4}p'-\frac{x}{2}(1-x^{2})^{-1/4}p\in L^{2}[-1,1],\\ (1-x^{2})^{-1/4}p \in L^{2}[-1,1]. $$ Equivalently, $$ \int_{-1}^{1}(1-x^{2})^{3/2}|p'|^{2}\,dx < \infty,\;\;\;\; \int_{-1}^{1}(1-x^{2})^{-1/2}|p|^{2}\,dx < \infty. $$ That means that the series solutions can diverge near the endpoints and still qualify as eigenfunctions, but not diverge too quickly. $p \sim (1-x^{2})^{-1/4}$ is borderline, and then $f=(1-x)^{1/2}p$ might be expected to vanish at $\pm 1$. That behavior should select out the correct eigenfunctions. One useful special case is $$ (1-x^{2})^{1/2}((1-x^{2})^{1/4}p)'\cdot (1-x^{2})^{-1/4}p \in L^{1} \\ \implies \frac{d}{dx}((1-x^{2})^{1/4}p)^{2}\in L^{1}\\ \implies (1-x^{2})^{1/2}p^{2} \mbox{ has limits at $x=\pm 1$ }. $$ These limits must be $0$ because $\frac{1}{1-x^{2}}((1-x^{2})^{1/2}p^{2})=((1-x^{2})^{-1/4}p)^{2}\in L^{1}$. So these conditions are necessary for the eigenfunctions $$ \lim_{x\rightarrow\pm 1}(1-x^{2})^{1/4}p(x)=0. $$