[proofplan]
We first isolate the Bernstein basis polynomials and observe that the defining formula is a finite linear combination of continuous functions, so $B_n f$ is again continuous. Linearity follows by applying the linearity of scalar multiplication and addition inside the finite sum. Positivity follows from the non-negativity of both the sampled values and the Bernstein basis functions on $[0,1]$. Finally, order preservation is reduced to positivity by applying $B_n$ to the non-negative function $g-f$.
[/proofplan]
[step:Define the Bernstein basis functions and check the operator is well-defined]
Fix $n \in \mathbb{N}$ with $n \geq 1$. For each integer $k$ with $0 \leq k \leq n$, define the Bernstein basis polynomial
\begin{align*}
b_{n,k}: [0,1] \to \mathbb{R}, \qquad b_{n,k}(x) = \binom{n}{k}x^k(1-x)^{n-k}.
\end{align*}
Each $b_{n,k}$ is the restriction to $[0,1]$ of a real polynomial, hence belongs to $C([0,1];\mathbb{R})$. For every $f \in C([0,1];\mathbb{R})$, the function $B_n f$ is the finite sum
\begin{align*}
B_n f = \sum_{k=0}^{n} f\left(\frac{k}{n}\right)b_{n,k}.
\end{align*}
Since each coefficient $f(k/n)$ is a real number and each $b_{n,k}$ is continuous on $[0,1]$, this finite sum belongs to $C([0,1];\mathbb{R})$. Therefore $B_n: C([0,1];\mathbb{R}) \to C([0,1];\mathbb{R})$ is well-defined.
[/step]
[step:Use the finite sum formula to prove linearity]
Let $f,g \in C([0,1];\mathbb{R})$ and let $\alpha,\beta \in \mathbb{R}$. For every $x \in [0,1]$, the definition of $B_n$ gives
\begin{align*}
(B_n(\alpha f+\beta g))(x) = \sum_{k=0}^{n}(\alpha f+\beta g)\left(\frac{k}{n}\right)b_{n,k}(x).
\end{align*}
Using the pointwise operations in $C([0,1];\mathbb{R})$, for each $k$ we have
\begin{align*}
(\alpha f+\beta g)\left(\frac{k}{n}\right)=\alpha f\left(\frac{k}{n}\right)+\beta g\left(\frac{k}{n}\right).
\end{align*}
Substituting this identity into the finite sum and distributing scalar multiplication over finite addition yields
\begin{align*}
(B_n(\alpha f+\beta g))(x)=\alpha (B_n f)(x)+\beta (B_n g)(x).
\end{align*}
Since this holds for every $x \in [0,1]$, we have
\begin{align*}
B_n(\alpha f+\beta g)=\alpha B_n f+\beta B_n g.
\end{align*}
Thus $B_n$ is linear.
[guided]
We prove linearity directly from the formula defining $B_n$. Let $f,g \in C([0,1];\mathbb{R})$ and let $\alpha,\beta \in \mathbb{R}$. The function $\alpha f+\beta g: [0,1]\to \mathbb{R}$ is defined pointwise by
\begin{align*}
(\alpha f+\beta g)(t)=\alpha f(t)+\beta g(t)
\end{align*}
for every $t \in [0,1]$. In particular, evaluating at the sampling point $k/n$ gives
\begin{align*}
(\alpha f+\beta g)\left(\frac{k}{n}\right)=\alpha f\left(\frac{k}{n}\right)+\beta g\left(\frac{k}{n}\right).
\end{align*}
Now fix $x \in [0,1]$. By the definition of the Bernstein operator and the basis functions $b_{n,k}$,
\begin{align*}
(B_n(\alpha f+\beta g))(x)=\sum_{k=0}^{n}(\alpha f+\beta g)\left(\frac{k}{n}\right)b_{n,k}(x).
\end{align*}
Replacing each sampled value by the pointwise formula above gives
\begin{align*}
(B_n(\alpha f+\beta g))(x)=\sum_{k=0}^{n}\left(\alpha f\left(\frac{k}{n}\right)+\beta g\left(\frac{k}{n}\right)\right)b_{n,k}(x).
\end{align*}
Because the sum is finite, we may distribute addition and scalar multiplication term by term:
\begin{align*}
(B_n(\alpha f+\beta g))(x)=\alpha\sum_{k=0}^{n}f\left(\frac{k}{n}\right)b_{n,k}(x)+\beta\sum_{k=0}^{n}g\left(\frac{k}{n}\right)b_{n,k}(x).
\end{align*}
The two sums on the right are precisely $(B_n f)(x)$ and $(B_n g)(x)$, so
\begin{align*}
(B_n(\alpha f+\beta g))(x)=\alpha (B_n f)(x)+\beta (B_n g)(x).
\end{align*}
Since the equality holds at every $x \in [0,1]$, the functions are equal. Therefore
\begin{align*}
B_n(\alpha f+\beta g)=\alpha B_n f+\beta B_n g.
\end{align*}
This is exactly the linearity of $B_n$.
[/guided]
[/step]
[step:Use non-negative coefficients to prove positivity]
Let $f \in C([0,1];\mathbb{R})$ satisfy $f(t)\geq 0$ for every $t \in [0,1]$. Fix $x \in [0,1]$. For each $k$ with $0\leq k\leq n$, the sampling point $k/n$ lies in $[0,1]$, so
\begin{align*}
f\left(\frac{k}{n}\right)\geq 0.
\end{align*}
Also $\binom{n}{k}\geq 0$, $x^k\geq 0$, and $(1-x)^{n-k}\geq 0$, hence
\begin{align*}
b_{n,k}(x)=\binom{n}{k}x^k(1-x)^{n-k}\geq 0.
\end{align*}
Therefore every term in the finite sum defining $(B_n f)(x)$ is non-negative:
\begin{align*}
f\left(\frac{k}{n}\right)b_{n,k}(x)\geq 0.
\end{align*}
A finite sum of non-negative [real numbers](/page/Real%20Numbers) is non-negative, so
\begin{align*}
(B_n f)(x)=\sum_{k=0}^{n}f\left(\frac{k}{n}\right)b_{n,k}(x)\geq 0.
\end{align*}
Since $x \in [0,1]$ was arbitrary, $B_n f\geq 0$ pointwise on $[0,1]$. Thus $B_n$ is positive.
[/step]
[step:Apply positivity to the difference to prove order preservation]
Let $f,g \in C([0,1];\mathbb{R})$ satisfy $f(x)\leq g(x)$ for every $x\in[0,1]$. Define
\begin{align*}
h: [0,1]\to\mathbb{R}, \qquad h(x)=g(x)-f(x).
\end{align*}
Since $f$ and $g$ are continuous, $h \in C([0,1];\mathbb{R})$. The assumed pointwise inequality gives $h(x)\geq 0$ for every $x\in[0,1]$. By positivity of $B_n$, proved above, $B_n h\geq 0$ pointwise on $[0,1]$. By linearity of $B_n$,
\begin{align*}
B_n h = B_n(g-f)=B_n g-B_n f.
\end{align*}
Hence, for every $x\in[0,1]$,
\begin{align*}
(B_n g)(x)-(B_n f)(x)\geq 0.
\end{align*}
Equivalently,
\begin{align*}
(B_n f)(x)\leq (B_n g)(x).
\end{align*}
Thus $B_n$ preserves pointwise order on $C([0,1];\mathbb{R})$, completing the proof.
[/step]