A famous idea is to approximate the integral $\int_{0}^{+\infty}e^{-x^2}\,dx$ with $\int_{0}^{+\infty}\frac{dt}{(1+t^2/n)^n}$ and $\int_{0}^{\sqrt{n}}(1-s^2/n)^n\,ds$, with $n\in\mathbb{N}$ tending to $+\infty$. These integrals are elementary: they can be computed through the substitutions $t=\sqrt{n}\tan\theta$, $s=\sqrt{n}\sin\varphi$ and repeated integration by parts.
The outcome is the double bound
$$ L(n)=\frac{\sqrt{n}4^n}{\binom{2n}{n}(2n+1)}\leq \int_{0}^{+\infty}e^{-x^2}\,dx \leq \frac{\pi n \sqrt{n}\binom{2n}{n}}{(2n-1) 4^n}=R(n)\tag{1}$$
holding for any $n\geq 1$. The statement
$$ \lim_{n\to +\infty}\frac{R(n)}{L(n)}=1\tag{2} $$
is equivalent to Wallis' product and the statement
$$ \lim_{n\to +\infty} R(n)L(n)=\frac{\pi}{4}\tag{3} $$
is trivial. By squeezing it follows that $\int_{0}^{+\infty}e^{-x^2}\,dx=\frac{\sqrt{\pi}}{2}$ and $\int_{-\infty}^{+\infty}e^{-x^2}\,dx = \sqrt{\pi}$.
Truth to be told, this is a bit of a fraud, too. We are not really avoiding the $\Gamma$ function, we are just relying on the reflection/duplication formulas for $\Gamma$ without making an explicit mention of $\Gamma$. Wallis' product itself is an instance of $\Gamma(s)\Gamma(1-s)=\frac{\pi}{\sin(\pi s)}$.
On the other hand, why to avoid the $\Gamma$ function? To say the least, it plays a major role in many relevant probability distributions, and as one of my menthors (C.Viola) likes to say, "a good mathematician or phycisist should not be afraid to manipulate the $\Gamma$ function or the sine function, also because they are not that different. The sooner one gets introduced to $\Gamma$, the better".
A (very!) shortened approach is to notice that $\Gamma\left(\frac{1}{2}\right)$ is the value we are looking for and $\Gamma\left(\frac{1}{2}\right)^2$ is the area of the unit circle.