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?
A weaker bound allows a completely elementary proof.
By Rodrigues formula we know that all the zeroes of $H_n$ are simple and real: let us denote them as $\zeta_1,\ldots,\zeta_n$. By Vieta's formulas the quantity $\sum_{k=1}^{n}\zeta_k^2$ only depends on the coefficients of $H_n$, and by the recurrence relations for Hermite polynomials we get
$$ \sum_{k=1}^{n}\zeta_k^2 =\frac{n(n-1)}{2},\qquad \sum_{k=1}^{n}\zeta_k^4 =\frac{n(n-1)(2n-3)}{4}, $$
$$ \sum_{k=1}^{n}\zeta_k^6 =\frac{n(n-1)(5n^2-17n+15)}{8},$$
$$\frac{7}{8}n(n-2)^4\leq\sum_{k=1}^{n}\zeta_k^8=\frac{n(n-1)(14n^3-79n^2+155n-105)}{16}\leq n\max_{k}|\zeta_k|^8 $$
$$\frac{7}{8}(n-1)^5\geq\sum_{k=1}^{n}\zeta_k^8=\frac{n(n-1)(14n^3-79n^2+155n-105)}{16}\geq 2\max_{k}|\zeta_k|^8 $$ so for any $n>2$ we have
$$\boxed{\frac{101}{112}(n-1)^{5/8}\geq \max_k |\zeta_k| \geq \frac{29}{30}\sqrt{n-2}.}$$
In order to prove the given upper bound one may use the integral representation for $H_n$ that comes from the residue theorem applied to the generating function: $$ H_n(x)=\frac{n!}{2\pi i}\oint_{\|z\|=\rho}\frac{e^{2xz-z^2}}{z^{n+1}}\,dz $$ Given some $x$ with a large absolute value ($\geq\sqrt{4n+1}$), it is enough to choose a suitable $\rho$ (close to $|x|$) and approximate the contour integral with decent accuracy, in order to show that its real part cannot be zero.
On the other hand, by quoting the answer of J.M. to the other question,
The high-minded, linear algebraic route involves deriving the recursion relation $$\hat{H}_{n+1}(x)=x\hat{H}_n(x)-\frac{n}{2}\hat{H}_{n-1}(x)$$ for the monic Hermite polynomial $\hat{H}_n(x)=2^{-n}H_n(x)$ (that is, the polynomial normalized to have unit leading coefficient), and from this derive the symmetric tridiagonal Jacobi matrix $$\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}$$ whose characteristic polynomial is $\hat{H}_n(x)$.
then applying the Gershgorin circle theorem, we immediately get $$ \max_k |\zeta_k| \leq \sqrt{2n-2} $$ for any $n>2$.
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}}. $$
