<s>[INST] Proving that all the real roots of Hermite polynomials are in $(-\sqrt{4n+1}, \sqrt{4n+1})$ The Hermite polynomials are given by:
$H_n(x)=(-1)^n e^{x^2} \dfrac{d^n}{dx^n}e^{-x^2}$
There is the proof that all the roots are real: https://math.stackexchange.com/a/104875/504137.
And I know the fact that they all are bounded i.e. all the roots lie in $(-\sqrt{4n+1}, \sqrt{4n+1})$. 
But how to prove that?
 [/INST] It doesn't look like there are many questions on this topic on Math.SE, so just for fun let's use some tricks from matrix analysis to slightly improve the bound $\sqrt{2n-2}$ in Jack's answer.

Let $A = [a_{ij}]$ be a real $n \times n$ matrix. Define $|A| = [|a_{ij}|]$. We will write $A \geq 0$ if all $a_{ij} \geq 0$. Let $\rho(A)$ be the absolute value of the largest eigenvalue of $A$.

Theorem. Let $A$ and $B$ be real $n \times n$ matrices. If $B - |A| \geq 0$ then $\rho(A) \leq \rho(|A|) \leq \rho(B)$.

This is theorem 8.1.18 in Horn & Johnson's Matrix Analysis.

Theorem. The eigenvalues of the $n \times n$ matrix
  $$
\begin{pmatrix}a&b&&&\\c&a&b&&\\&c&\ddots&\ddots&\\&&\ddots&\ddots&b\\&&&c&a\end{pmatrix}
$$
  are
  $$
\lambda_k = a + 2\sqrt{bc} \cos\left(\frac{k\pi}{n+1}\right), \qquad k=1,\ldots,n.
$$

This is a known fact about tridiagonal Toeplitz matrices. See, e.g., this PDF.

Using J.M.'s matrix from Jack's answer, the matrix
$$
\begin{pmatrix}0&\sqrt{\frac{n-1}2}&&&\\\sqrt{\frac{n-1}2}&0&\sqrt{\frac{n-1}2}&&\\&\sqrt{\frac{n-1}2}&\ddots&\ddots&\\&&\ddots&\ddots&\sqrt{\frac{n-1}2}\\&&&\sqrt{\frac{n-1}2}&0\end{pmatrix} -
 \begin{pmatrix}0&\sqrt{\frac12}&&&\\\sqrt{\frac12}&0&\sqrt{\frac22}&&\\&\sqrt{\frac22}&\ddots&\ddots&\\&&\ddots&\ddots&\sqrt{\frac{n-1}2}\\&&&\sqrt{\frac{n-1}2}&0\end{pmatrix}
$$
is $\geq 0$ for $n \geq 2$, so, by the two theorems above,
$$
\max_k |\zeta_k| \leq \sqrt{2n-2} \cos\left(\frac{\pi}{n+1}\right).
$$
To get a lower bound we can note that when $n$ is even the characteristic polynomial of the matrix
$$
\begin{pmatrix}
0&\sqrt{\frac12}\\
\sqrt{\frac12}&0&0\\
&0&0&\sqrt\frac32\\
&&\sqrt\frac32&0&0\\
&&&0&0&\sqrt\frac52\\
&&&&\sqrt\frac52&0&\ddots\\
&&&&&\ddots&\ddots&\sqrt\frac{n-1}2\\
&&&&&&\sqrt\frac{n-1}2&0
\end{pmatrix}
$$
is
$$
\prod_{k=0}^{(n-2)/2} \left(\lambda^2 - \frac{2k+1}{2}\right),
$$
and when $n$ is odd the characteristic polynomial of the matrix
$$
\begin{pmatrix}
0&0\\
0&0&\sqrt\frac22\\
&\sqrt\frac22&0&0\\
&&0&0&\sqrt\frac42\\
&&&\sqrt\frac42&0&0\\
&&&&0&0&\ddots\\
&&&&&\ddots&\ddots&\sqrt\frac{n-1}2\\
&&&&&&\sqrt\frac{n-1}2&0
\end{pmatrix}
$$
is
$$
-\lambda \prod_{k=0}^{(n-1)/2} \left(\lambda^2 - k\right).
$$
Since these two matrices are $\leq$ the Hermite matrix, we get the lower bound
$$
\max_k |\zeta_k| \geq \sqrt{\frac{n-1}{2}}.
$$
</s>