How to show simply that $e^{\frac{x}{2}}\int^\infty_0 e^{-t}t^{n-\frac{1}{2}}\cos(2\sqrt{xt})dt=O\big(\frac{n!}{\sqrt{n}}\big)$?

Making the normalised change of variables $t = ns^2$, $x = 4ny^2$ (with $y \geq 0$ and $s$ of either sign) one can write $$ f_n(x) = e^{2ny^2} n^n \sqrt{n} \int_{{\bf R}} e^{-n\phi(s)}\ ds$$ where the phase $\phi(s)$ is given by $$ \phi(s) := s^2 - 2 \log s - 4iy s$$ using the standard branch of the complex logarithm. This phase has stationary points at $iy \pm \sqrt{1-y^2}$. In the "Bessel" regime $y < 1-\delta$ for some fixed $\delta>0$, the stationary points are non-degenerate, and $\mathrm{Re} \phi(s)$ attains a local minimum of $1 + 2y^2$ at these points (with respect to a horizontal contour), so one gets the required bound $f_n(x) \lesssim n^n/e^n$ (equivalent to the claimed bound by Stirling's formula) in this case by the usual saddle point method. Similarly, in the "exponential" regime $y > 1 + \delta$, the stationary points are again non-degenerate; if one makes the convenient substitution $y = \cosh \theta$ with $\theta>0$ then at the stationary point $i e^{\theta}$, $\mathrm{Re} \phi(s)$ attains a local minimum (again wrt a horizontal contour) of $1 + 2y^2 + (\sinh 2\theta - 2\theta) > 1 + 2y^2$, and stationary phase again gives the desired bound (with room to spare). However in the "Airy" regime in which $y$ is close to $1$ (or in your original coordinates, $x$ is close to $4n$) the stationary points coalesce and I don't think your claimed bounds actually hold (one can lose an additional factor of $n^{1/6}$, I think, and this is consistent with the standard Airy asymptotics given for instance here). You may want to double check your reference for the Laguerre polynomial bounds; they may not hold in the Airy regime $x \approx 4n$. For instance, the bound in Theorem 8.91.2 of Szego's book only gives a decay bound of $n^{-1/3}$ instead of $n^{-1/2}$ once one works in a range of parameters that includes the Airy regime.