This method might not be efficient enough. Denote $ f_n(x)$ as the function with degree $ n$, we can have $ f_n(x)=g_n(x)+h_n(x)$, where $ g_n$ is even and $ h_n$ is odd. $ (x+1)f(x)^2-1$ being odd implies $ (1+x)f(x)^2+(1-x)f(-x)^2=2\Rightarrow g(x)^2+h(x)^2+2xg(x)h(x)=1$ (*) Now we define $ g_n,h_n$ as follows: $ g_0(x)=1$,$ h_0(x)=0$. $ g_{2n+1}(x)=g_{2n}(x)$, $ h_{2n+1}(x)=-2xg_{2n}(x)-h_{2n}(x)$; $ g_{2n}(x)=g_{2n-1}(x)+2xh_{2n-1}(x)$, $ h_{2n}(x)=-h_{2n-1}(x)$. It is not hard to verify (*), $ g_n(0)=1$, $ g_n$ being even, and $ h_n$ being odd. We proved existence. Now we prove uniqueness by induction. Uniqueness of $ f_0(x)$ is trivial. Suppose $ f_{n-1}(x)$ is unique and we look at $ f_n(x)=g_n(x)+h_n(x)$. Assume $ n$ is even. $ (g_n,h_n)$ is any pair which meets the requirement. Then we see $ g_n$ has degree $ n$. We also have $ g_n^2+h_n^2+2xg_nh_n=1$. If we look at the coefficient of $ x^{2n}$ of the above equation, we can show $ h_n$ has degree $ n-1$ and $ g_n+2xh_n$ has degree no greater than $ n-2$. Now we find out $ g=g_n+2xh_n$, $ h=-h_n$ also solves (*). But we notice $ g+h$ has degree $ n-1$. Then by induction assumption, we have $ g_n+2xh_n=g_{n-1}$ and $ -h_n=h_{n-1}$. This shows uniqueness of $ f_n$ Similarly we can handle the case when $ n$ is odd. Q.E.D.

All I want to say is:a perfect solution! Siince the official solution is based on Pell Equation(defined in polynomial not integer) which is not elementary at all

Could you post the official solution of this problem,dear hxy09?

We will show that up to sign, there is in fact a unique polynomial in $\mathbb{R}[x]$. Let $f(x)=p(x^2)+xq(x^2)$, so $(x+1)f(x)^2-1=-[(-x+1)f(-x)^2-1]$ rewrites as $x[p(x)+q(x)]^2+(1-x)p(x)^2=1$. (*) Let $r(x)=p(x)+q(x)$. Differentiating both sides, we have \[r(x)[r(x)+2xr'(x)]=p(x)[p(x)+2(x-1)p'(x)].\]Clearly $\deg{p}=\deg{r}$ or else (*) can't hold (the leading coefficients must be the same up to sign), so because $\gcd(p(x),r(x))=1$, $p(x)$ must divide $r(x)+2xr'(x)$ and so there exists a constant $c$ such that $p(x)=c[r(x)+2xr'(x)]$. Plugging this into (*), we get the differential equation \[1=[x+c^2(1-x)]r(x)^2+4c^2x(1-x)r'(x)r(x)+4c^2x^2(1-x)r'(x)^2.\]From the quadratic formula, we find \[r'(x)=-\frac{c(x-1)r(x)\pm\sqrt{(x-1)[xr(x)^2-1]}}{2cx(x-1)},\]or letting $s(x)=r(x)\sqrt{x}$ so that $2xr'(x)\sqrt{x}=2xs'(x)-s(x)$ and simplifying, \[\frac{2cs'(x)}{\sqrt{s(x)^2-1}}=\sqrt{\frac{1}{(2x-1)^2-1}},\]where WLOG $\lim_{x\to\infty}s(x)=\infty$. Anyway, it's easy to find that \[\int\frac{1}{\sqrt{x^2-1}}\;dx=\ln(x+\sqrt{x^2-1})+C,\]so for some constant $C_0$ we have \[c\ln(s(x)+\sqrt{s(x)^2-1})=\ln(\sqrt{x}+\sqrt{x-1})+C_0,\]i.e. \[\frac{1}{s(x)-\sqrt{s(x)^2-1}}=s(x)+\sqrt{s(x)^2-1}=e^{C_0/c}(\sqrt{x}+\sqrt{x-1})^{1/c},\]so there exist constants $C_1,C_2$ such that \[r(x)=\frac{e^{C_1}(\sqrt{x}+\sqrt{x-1})^{C_2}+e^{-C_1}(\sqrt{x}-\sqrt{x-1})^{C_2}}{2\sqrt{x}}.\]Since $r$ is a polynomial, we clearly have $C_2\in\mathbb{Z}$ and thus $C_1=0$. Now we can easily compute $f(x)=(1-x)p(x^2)+xr(x^2)=T_{n+1}(x)-(x-1)U_n(x)$ (up to sign of course). Note that (*) rearranges to the Pell equation $u(x)^2-(x^2-1)v(x)^2=1$ upon the substitutions $u(x)=xr(x^2)$, $v(x)=p(x^2)$, which is a bit cleaner to do with differentiation.

Solution from Twitch Solves ISL: Suppose the polynomial $P$ obeys \[ 2 = (x+1)P(x)^2 + (-x+1)P(-x)^2. \]Work in ${\mathbb R}[y, y^{-1}]$ and substitute \[ x = \frac{1}{2}\left( y + \frac 1y \right). \]Then the above equation can be rewritten as \begin{align*} 4 &= \left( y + \frac 1y + 2 \right) P(x)^2 - \left( y + \frac 1y - 2 \right) P(-x)^2 \\ \iff 4y &= \left( (y+1) P(x) \right)^2 - \left( (y-1) P(-x) \right)^2 \\ \iff 4y &= \Big( (y+1) P(x) + (y-1) P(-x) \Big) \Big( (y+1) P(x) - (y-1) P(-x) \Big). \end{align*}Using factorization in ${\mathbb R}[y, y^{-1}]$, we see there must exist an integer $d \ge 0$ and a nonzero real number $\varepsilon$ such that we have identities \begin{align*} 2\varepsilon \cdot y^{d+1} &= (y+1) P(x) \pm (y-1) P(-x) \\ \frac 2\varepsilon \cdot y^{-d} &= (y+1) P(x) \mp (y-1) P(-x) \end{align*}for some opposite choice of signs $\pm$. Replacing $y$ with $1/y$, it follows $\varepsilon = \varepsilon^{-1} \implies \varepsilon = \pm 1$. So the system of equations is equivalent to \begin{align*} y^{d+1} + y^{-d} &= \varepsilon (y+1) P(x) \\ y^{d+1} - y^{-d} &= \pm \varepsilon (y-1) P(-x) \end{align*}Now, the first equation implies that \begin{align*} \varepsilon P(x) &= \frac{y^{d+1}+y^{-d}}{y+1} \\ &= \left( y^d + \frac{1}{y^d} \right) - \left( y^{d-1} + \frac{1}{y^{d-1}} \right) + \dots + (-1)^{d-1} \left( y + \frac 1y \right) + (-1)^d. \end{align*}The theory of Chebyshev polynomials guarantees that $y^k + \frac{1}{y^k}$ is a degree-$k$ polynomial in $x$ for each $k \ge 0$. Hence we see there is a unique choice of polynomial $P$ of degree $n$ (corresponding to $d=n$) which makes the identity $P(x) = \varepsilon \frac{y^{d+1}+y^{-d}}{y+1}$ true, for each $\varepsilon \in \{\pm1\}$. When we add the condition $P(0) = 1$ (set $y=i$ so $x=0$), we find only one $\varepsilon \in \{\pm1\}$ is valid. Conversely, we have the implication \[ \varepsilon (y+1)P(x) = y^{d+1} + y^{-d} \implies \varepsilon (y-1)P(-x) = \pm (y^{d+1} - y^{-d}) \]which follows readily by making the substitution $y \mapsto iy$. So the choice of polynomial $P$ we just mentioned will satisfy the entire system, ergo (reversing the logic above) will satisfy the desired condition. In other words, the unique possibility for $P$ does indeed work. Remark: The first few actual polynomials $P$ are: \begin{align*} P(x) &= 1 \\ P(x) &= 1 - 2x \\ P(x) &= 1 + 2x - 4x^2 \\ P(x) &= 1 - 4x - 4x^2 + 8x^3 \\ P(x) &= 1 + 4x - 12x^2 - 8x^3 + 16x^4. \end{align*}