Browsing through my Zotero database, I fall upon a paper by Harald Kuhn where he proposes an elementary proof of the Weierstrass approximation theorem. The proof is indeed neat, so here it is.
Theorem. — Let $f\colon[0;1]\to\mathbf R$ be a continuous function and let $\varepsilon$ be a strictly positive real number. There exists a polynomial $P\in\mathbf R[T]$ such that $\left|P(x)-f(x)\right|<\varepsilon$ for every $x \in[0;1]$.
The proof goes out of the world of continuous functions and considers the Heaviside function $H$ on $\mathbf R$ defined by $H(x)=0$ for $x<0$ and $H(x)=1$ for $x\geq 0$ and will construct a polynomial approximation of $H.$
Lemma. — There exists a sequence $(P_n)$ of polynomials which are increasing on $[-1;1]$ and which, for every $\delta>0$, converge uniformly to the function $H$ on $[-1;1]\setminus [-\delta;\delta]$.
For $n\in\mathbf N$, consider the polynomial $Q_n=(1-T^n)^{2^n}$. On $[0;1]$, it defines an decreasing function, with $Q_n(0)=1$ and $Q_n(1)=0$.
Let $q\in[0;1]$ be such that $0\leq q<1/2$. One has $$ Q_n(q) = (1-q^n)^{2^n} \geq 1-2^n q^n $$ in view of the inequality $(1-t)^n \geq 1-nt$ which is valid for $t\in[0;1]$. Since $2q<1$, we see that $Q_n(q)\to 1$. Since $Q_n$ is decreasing, one has $Q_n(q)\leq Q_n(x)\leq Q_n(1)=1$ for every $x\in[0;q]$, which shows that $Q_n$ converges uniformly to the constant function $1$ on the interval $[0;q]$.
Let now $q\in[0;1]$ be such that $1/2<q\leq 1$. Then $$ \frac1{Q_n(q)} = \frac1{(1-q^n)^{2^n}} = \left(1 + \frac{q^n}{1-q^{n}}\right)^{2^n} \geq 1 + \frac{2^nq^n}{1-q^{n}} $$ so that $Q_n(q)\to 0$. Since $Q_n$ is decreasing, one has $0=Q_n(1)\leq Q_n(x)\leq Q_n(q)$ for every $x\in[q;1]$, so that $Q_n$ converges uniformly to the constant function $0$ on the interval $[q;1]$.
Make an affine change of variables and set $P_n = Q_n((1-T)/2)$. The function defined by $P_n$ on $[-1:1]$ is increasing; for any real number $\delta$ such that $\delta>0,$ it converges uniformly to $0$ on $[-1;-\delta]$, and it converges uniformly to $1$ on $[\delta;1]$. This concludes the proof of the lemma.
We can now turn to the proof of the Weierstrass approximation theorem. Let $f$ be a continuous function on $[0;1]$. We may assume that $f(0)=0.$
The first step, that follows from Heine's uniform continuity theorem, consists in noting that there exists a uniform approximation of $f$ by a function of the form $ F(x)=\sum_{m=1}^N a_m H(x-c_m)$, where $(a_1,\dots,a_m)$ and $(c_1,\dots,c_m)$ are real numbers in $[0;1].$ Namely, Heine's theorem implies that there exists $\delta>0$ such that $|f(x)-f(y)|<\varepsilon$ if $|x-y|<\delta$. Choose $N$ such that $N\delta\geq 1$ and set $c_m=m/N$; then define $a_1,\dots,a_N$ so that $a_1=f(c_1)$, $a_1+a_2=f(c_2)$, etc. It is easy to check that $|F(x)-f(x)|\leq \varepsilon$ for every $x\in[0;1]$. Moreover, $\lvert a_m\rvert\leq\varepsilon/2$ for all $m$.
Now, fix a real number $\delta>0$ small enough so that the intervals $[c_m-\delta;c_m+\delta]$ are pairwise disjoint, and $n\in\mathbf N$ large enough so that $|P_n(x)-H(x)|\leq\varepsilon/2A$ for all $x\in[-1;1]$ such that $\delta\leq|x|$, where $A=\lvert a_0\rvert+\dots+\lvert a_N\rvert$. Finally, set $P(T) = \sum_{m=1}^N a_m P_n(T-c_m)$.
Let $x\in[-1;1]$. If $x$ doesn't belong to any interval of the form $[c_m-\delta;c_m+\delta]$, one can write $$\lvert P(x)-F(x)\rvert\leq \sum_{m} \lvert a_m\rvert \,\lvert P_n(x-c_m)- H(x-c_m)\rvert \leq \sum_m \lvert a_m\rvert (\varepsilon/A)\leq \varepsilon. $$ On the other hand, if there exists $m\in\{1,\dots,N\}$ such that $x\in[c_m-\delta;c_m+\delta]$, then there exists a unique such integer $m$. Writing $$\lvert P(x)-F(x)\rvert\leq \sum_{k\neq m} \lvert a_k\rvert \,\lvert P_n(x-c_k)- H(x-c_k)\rvert + \lvert a_m\rvert\, \lvert P_n(x-c_m)-H(x-c_m)\rvert, $$ the term with index $k$ in the first sum is bounded by $\lvert a_k\rvert \varepsilon/A$, while the last term is bounded by $ \lvert a_m\rvert$, because $0\leq P_n\leq H\leq 1$. Consequently, $$\lvert P(x)-F(x)\rvert \leq (\varepsilon/2A)\sum_{k\neq m} \lvert a_k\rvert +\lvert a_m\rvert \leq 2\varepsilon.$$ Finally, $\lvert P(x)-f(x)\rvert\leq \lvert P(x)-F(x)\rvert+\lvert F(x)-f(x)\rvert\leq 3\varepsilon$. This concludes the proof.