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?

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}}.
$$
