The following proof has two parts: in the first part I will prove that $n!\sim C\sqrt{n}\left(\frac{n}{e}\right)^n$, in the second part I will show that $C=\sqrt{2\pi}$. It is crucial to notice that
$$ m = \prod_{k=1}^{m-1}\left(1+\frac{1}{k}\right)\tag{1}$$
from which:
$$ n! = \prod_{m=2}^{n}\prod_{k=1}^{m-1}\left(1+\frac{1}{k}\right) = \prod_{k=1}^{n-1}\left(1+\frac{1}{k}\right)^{n-k}=\frac{n^n}{\prod_{k=1}^{n-1}\left(1+\frac{1}{k}\right)^k}\tag{2}$$
and
$$ \log\left(\frac{n^n}{n!}\right) = (n-1)-\sum_{k=1}^{n-1}\left[1-k\log\left(1+\frac{1}{k}\right)\right] \tag{3} $$
where $1-k\log\left(1+\frac{1}{k}\right)$ is bounded between $\frac{1}{2k}-\frac{1}{3k^2}$ and $\frac{1}{2k}$ for any $k\geq 1$. This leads to
$$ \log\left(\frac{n^n}{n!}\right) = n-\frac{1}{2}\log(n)+K+O\left(\frac{1}{n}\right)\tag{4} $$
and the first part is granted, $n!\sim C\sqrt{n}\left(\frac{n}{e}\right)^n$. As a consequence,
$$ \lim_{n\to +\infty}\frac{\sqrt{n}}{4^n}\binom{2n}{n} = \lim_{n\to +\infty}\frac{C\sqrt{2n}\left(\frac{2n}{e}\right)^{2n}\sqrt{n}}{4^n C^2 n \left(\frac{n}{e}\right)^{2n}}=\frac{\sqrt{2}}{C}\tag{5}$$
and it is enough to show that $ \lim_{n\to +\infty}\frac{\sqrt{n}}{4^n}\binom{2n}{n} = \frac{1}{\sqrt{\pi}}$. On the other hand
$$ \frac{(2n+1)!!}{(2n)!!}=\prod_{k=1}^{n}\left(1+\frac{1}{2k}\right) = \sqrt{\prod_{k=1}^{n}\left(1+\frac{1}{k}\right)\prod_{k=1}^{n}\left(1-\frac{1}{(2k+1)^2}\right)^{-1}} \tag{6} $$
leads to the identity
$$ \frac{2n+1}{4^n}\binom{2n}{n} = \sqrt{n+1}\sqrt{\prod_{k=1}^{n}\left(1-\frac{1}{(2k+1)^2}\right)^{-1}} \tag{7} $$
and
$$ \prod_{k\geq 1}\left(1-\frac{1}{(2k+1)^2}\right)=\frac{\pi}{4}\tag{8} $$
by the Weierstrass product for the cosine function. This finishes the proof.
As an alternative,
$$ \frac{(2n+1)!!}{(2n)!!}=\prod_{k=1}^{n}\left(1+\frac{1}{k}\right) = \frac{2\,\Gamma\left(\frac{3}{2}+n\right)}{\Gamma\left(\frac{1}{2}\right)\,\Gamma(n+1)}\sim \frac{2}{\Gamma\left(\frac{1}{2}\right)}\sqrt{n}\tag{9} $$
for large $n$s by Gautschi's inequality ($\log\Gamma$ is a convex function on $\mathbb{R}^+$, by the Bohr-Mollerup characterization of just by applying the Cauchy-Schwarz inequality to its integral representation), and
$$ \Gamma\left(\frac{1}{2}\right) = \int_{0}^{+\infty}z^{-1/2}e^{-z}\,dz = 2 \int_{0}^{+\infty}e^{-x^2}\,dx = \int_{-\infty}^{+\infty}e^{-x^2}\,dx = \sqrt{\pi}.\tag{10}$$