That's nice! Thanks for your careful answer. :)

And I have a question for this answer, how to guarantee $\sqrt{x(x+4)}=t^2-\frac{1}{t^2}$ rather than $\sqrt{x(x+4)}=-t^2+\frac{1}{t^2}$, and so on.

Thanks, but the interval of the root of $f_n$ seems to be $\left [ -4,0 \right ]$?

@Per Alexandersson I tried to $n=200$, it still holds.

Thank you for your answer. We have also tried this method, but the sign pattern of $P_n$ is $--++\cdots $, and the existing results cannot seem to get $P_n+Q_n \in RZ$

Thanks for Timothy Chow and Ira Gessel's answers, amazing construction! At first, we tried to use the Sturm sequence to prove the real zeros, but it seem to be impossible :(

But in this recurrence relation, we can't seem to be sure that $f_n(x)$ is all real roots.

Thanks for your comment! I'm sorry that I forgot to indicate the $n \geq 1$. Based on your comment, then by the generating function of $f_n(x)$, we get the recurrence relation of $f_n(x)$, that is $f_n(x)=xf_{n-1}(x)+x^2f_{n-3}+x^2f_{n-4}$ (for $n \geq 5$).

Thanks for your answer, but I want to continue solving the determinant composed of eigenvectors, which seems difficult.