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.
