The needed observation here is that, whereas the relative error of the Stirling approximation is worse for large $|x|$, it gets multiplied by a fast decreasing normal density function, and then the product gets integrated to produce the desired result. Here are the details.
Let $X$ be a binomial random variable with parameters $n\in\mathbb{N}$ and $p\in(0,1)$, where $n\to\infty$ and $p$ is fixed. For any integer $k$, let
$$z_k:=\frac{k-np}{\sqrt{npq}},\quad q:=1-p.
$$
Let
$$K:=\{k=0,\dots,n\colon|z_k|\le n^{1/4}\},$$
so that $k/n\to p$ uniformly over $k\in K$.
Let
$$\hat p:=k/n,\quad\hat q:=1-\hat p=(n-k)/n,$$
$$\epsilon:=\epsilon_k:=z_k/\sqrt n,\quad\lambda:=\sqrt{q/p},
$$
whence
$$\frac{\hat p}{p}=1+\epsilon\lambda,\quad\frac{\hat q}{q}=1-\epsilon/\lambda.$$
Uniformly over $k\in K$ one has $k\sim np$ and $n-k\sim nq$, so that $\frac1k=O(\frac1n)$ and $\frac1{n-k}=O(\frac1n)$.
Here and in what follows, the constant in $O(\cdot)$ depends only on $p$. So, by Stirling's formula, over all $k\in K$
$$P(X=k)=\frac{n!}{k!(n-k)!} p^kq^{n-k}
=\frac1{\sqrt{2\pi n\hat{p}\hat{q}}}\,\frac1R\,\Big(1+O\Big(\frac1n\Big)\Big)$$
$$
=\frac1{\sqrt{2\pi npq}}\,\frac1R\,\Big(1+O\Big(|\epsilon|+\frac1n\Big)\Big)
=\frac\delta{\sqrt{2\pi}}\,e^{-nH}\Big(1+O\Big(|\epsilon|+\frac1n\Big)\Big)
$$
where
$$R:=\Big(\frac{\hat p}{p}\Big)^k \Big(\frac{\hat q}{q}\Big)^{n-k},
$$
$$\delta:=\frac1{\sqrt{npq}},\quad H:=H_k:=\ln(R^{1/n})=\hat p\ln\frac{\hat p}p+\hat q\ln\frac{\hat q}q
:=p\psi(\epsilon\lambda)+q\psi(-\epsilon/\lambda),$$
$$\psi(s):=(1+s)\ln(1+s).
$$
Since $\psi(s)=s + s^2/2 +O(|s|^3)$ for $|s|\le1/2$ and
$\frac{|z_k|^3}{\sqrt n}+|\epsilon|=\frac{|z_k|^3}{\sqrt n}+\frac{|z_k|}{\sqrt n}=O\big(\frac{|z_k|^3+1}{\sqrt n}\big)$, we find that
$$(1)\qquad P(X=k)=\delta\,\phi(z_k)\Big(1+O\Big(\frac{|z_k|^3}{\sqrt n}\Big)\Big)\Big(1+O\Big(|\epsilon|+\frac1n\Big)\Big)$$
$$
=\delta\,\phi(z_k)\Big(1+O\Big(\frac{|z_k|^3+1}{\sqrt n}\Big)\Big)
$$
over all $k\in K$,
where $\phi$ is the standard normal density.
Next,
for $u=O(n^{1/4})$ and $|h|\le\delta/2$,
$$\frac{\phi(u+h)+\phi(u-h)}{2}=\phi(u)\cosh(hu)e^{-h^2/2}
=\phi(u)(1+O(h^2u^2+h^2))=\phi(u)\Big(1+O\Big(\frac1{\sqrt n}\Big)\Big).
$$
So, for $k\in K$
$$\delta\,\phi(z_k)=\frac12\,\int_{-\delta/2}^{\delta/2}[\phi(z_k+h)+\phi(z_k-h)]\,dh\,\Big(1+O\Big(\frac1{\sqrt n}\Big)\Big)$$
$$
=\int_{z_k-\delta/2}^{z_k+\delta/2}\phi(z)\,dz\,\Big(1+O\Big(\frac1{\sqrt n}\Big)\Big).
$$
Also,
for any $k\in K$ and any $z\in[z_k-\delta/2,z_k+\delta/2]$ one has
$|z_k|^3-|z|^3\le3\delta(|z_k|+\delta/2)^2=O(1)$, whence $|z_k|^3/\sqrt n=O((|z|^3+1)/\sqrt n)$.
So, by $(1)$, for some positive real $c$ depending only on $p$ and for all $k\in K$
$${J_k\!}^-\le P(X=k)\le J_k^+,$$
where
$$J_k^\pm:=\int_{z_k-\delta/2}^{z_k+\delta/2}\phi(z)
\Big(1\pm\frac{c(|z|^3+1)}{\sqrt n}\Big)\,dz.
$$
Hence, for any $x\in[0,n^{1/4}]$, letting $K_x:=\{k\in\mathbb{Z}\colon0<z_k\le x\}$, one has
$$P(0<B_n\le x)=\sum_{k\in K_x}P(X=k)
\le\sum_{k\in K_x}J_k^+$$
$$
\le\int_{-\delta/2}^{x+\delta/2}\phi(z)\Big(1+\frac{c(|z|^3+1)}{\sqrt n}\Big)\,dz$$
$$
\le\int_0^x\phi(z)\,dz+A\delta+\frac B{\sqrt n}
\le\int_0^x\phi(z)\,dz+\frac C{\sqrt n},
$$
where $A,B,C$ depend only on $p$. Quite similarly one obtains the corresponding lower bound on $P(0<B_n\le x)$, so that
$$(2)\qquad |P(0<B_n\le x)-(\Phi(x)-1/2)|\le\frac C{\sqrt n}
$$
for $x\in[0,n^{1/4}]$, where $\Phi$ is the standard normal cumulative distribution function.
On the other hand, for real $x>n^{1/4}$ one has $1-\Phi(x)=O(1/\sqrt n)$ and,
by Bernstein's inequality,
$$P(B_n>x)\le\exp\Big(-\frac{x^2}{2+x/\sqrt{npq}}\Big)
\le e^{-x^2/4}+e^{-x\sqrt{npq}/2}=O(1/\sqrt n).
$$
It follows that $(2)$ holds for all real $x\ge0$ (perhaps with a greater constant $C$). Similarly, $|P(B_n\le0)-1/2)|\le\frac C{\sqrt n}$.
Thus, for all real $x$
$$|P(B_n\le x)-\Phi(x)|\le\frac{2C}{\sqrt n},
$$
as desired.
Best Answer
Use concentration of measure. For 2 dimensions, Hoeffding's inequality. See the appendix of the book on empirical processes by van der Vaart and Wellenr for the multinomial case.