[proofplan]
We construct $d_A$ locally on each trivialising chart by the formula $d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha$, where $\theta_\alpha$ is the connection matrix in that chart, and then verify that the local definitions are compatible on overlaps. Compatibility reduces to a direct computation using the transformation laws for the local representative $\omega_\alpha$ and the connection matrix $\theta_\alpha$ under change of trivialisation. Once local consistency is established, the globally defined $E$-valued $(i+1)$-form $d_A \omega$ exists by the sheaf property of sections, and $\mathbb{R}$-linearity of the assignment $\omega \mapsto d_A \omega$ follows from linearity of each local formula. For $i = 0$ the operator restricts to the original connection, so this construction extends $d_A$ consistently to all degrees.
[/proofplan]
[step:Fix a trivialising atlas and recall the transition laws for local representatives and connection matrices]
Let $\{(U_\alpha, \psi_\alpha)\}_{\alpha \in I}$ be a trivialising atlas for $E$ with transition functions $\psi_{\alpha\beta}: U_\alpha \cap U_\beta \to \operatorname{GL}_r(\mathbb{R})$, and let $\theta_\alpha \in \Omega^1(U_\alpha; \mathrm{End}(\mathbb{R}^r))$ be the connection matrices of $A$ in each chart. Recall the transition laws:
1. For $\omega \in \Omega^i(E)$ with local representatives $\omega_\alpha \in \Omega^i(U_\alpha; \mathbb{R}^r)$ (obtained by applying the frame isomorphism $\psi_\alpha$), the overlap compatibility is
\begin{align*}
\omega_\alpha = \psi_{\alpha\beta}\, \omega_\beta \qquad \text{on } U_\alpha \cap U_\beta.
\end{align*}
2. The connection matrices transform as
\begin{align*}
\theta_\alpha = \psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} + \psi_{\alpha\beta}\, d(\psi_{\alpha\beta}^{-1}) = \psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} - d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1},
\end{align*}
using $d(\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}) = 0$ to rewrite $\psi_{\alpha\beta}\, d(\psi_{\alpha\beta}^{-1}) = -d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}$.
[/step]
[step:Define the local candidate $d_A \omega$ on each chart by $d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha$]
For $\omega \in \Omega^i(E)$, define locally
\begin{align*}
(d_A \omega)_\alpha := d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha \;\in\; \Omega^{i+1}(U_\alpha; \mathbb{R}^r),
\end{align*}
where $d$ denotes the component-wise de Rham differential on $\mathbb{R}^r$-valued forms and $\theta_\alpha \wedge \omega_\alpha$ is the wedge product combined with matrix--vector multiplication: $(\theta_\alpha \wedge \omega_\alpha)_j = \sum_{k} (\theta_\alpha)_{jk} \wedge (\omega_\alpha)_k$. For $i = 0$ this reduces to the connection formula $d_A s = ds + \theta_\alpha \cdot s$.
[guided]
The formula $d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha$ generalises the connection formula $d_A s = ds + \theta_\alpha \cdot s$ from sections (degree $i = 0$) to $E$-valued forms of arbitrary degree.
Let us unpack the notation. An element $\omega \in \Omega^i(E) = \Gamma(E \otimes \Lambda^i T^*M)$ has, in the trivialisation over $U_\alpha$, a local representative $\omega_\alpha \in \Omega^i(U_\alpha; \mathbb{R}^r)$: it is a vector of $i$-forms, one for each component of the local frame. The de Rham differential $d$ acts component-wise: $(d\omega_\alpha)_j = d((\omega_\alpha)_j)$ for $j = 1, \ldots, r$, and produces an $(i+1)$-form-valued vector in $\Omega^{i+1}(U_\alpha; \mathbb{R}^r)$.
The second term $\theta_\alpha \wedge \omega_\alpha$ is the wedge combined with matrix multiplication: $\theta_\alpha$ is a matrix of 1-forms, $\omega_\alpha$ is a vector of $i$-forms, and
\begin{align*}
(\theta_\alpha \wedge \omega_\alpha)_j := \sum_{k=1}^r (\theta_\alpha)_{jk} \wedge (\omega_\alpha)_k \in \Omega^{i+1}(U_\alpha),
\end{align*}
producing another $(i+1)$-form-valued vector. So both terms live in $\Omega^{i+1}(U_\alpha; \mathbb{R}^r)$ and their sum is well-defined.
Why this particular formula? For $i = 0$ sections it is literally the definition of the connection $d_A$ in the frame. For higher $i$ it is determined by the requirement that $d_A$ satisfy a graded Leibniz rule $d_A(\omega \wedge \eta) = d\omega \wedge \eta + (-1)^i \omega \wedge d_A \eta$ when $\omega$ is a scalar form, which fixes it uniquely once we know it on sections.
[/guided]
[/step]
[step:Verify compatibility on overlaps by computing $(d_A \omega)_\alpha - \psi_{\alpha\beta} (d_A \omega)_\beta$]
We claim that on $U_\alpha \cap U_\beta$,
\begin{align*}
(d_A \omega)_\alpha = \psi_{\alpha\beta}\, (d_A \omega)_\beta.
\end{align*}
Start from $\omega_\alpha = \psi_{\alpha\beta}\, \omega_\beta$ and apply $d$:
\begin{align*}
d\omega_\alpha = d(\psi_{\alpha\beta}\, \omega_\beta) = d\psi_{\alpha\beta} \wedge \omega_\beta + \psi_{\alpha\beta}\, d\omega_\beta,
\end{align*}
using the matrix Leibniz rule (the sign is positive because $\psi_{\alpha\beta}$ is a 0-form).
For the second term, using the transformation law $\theta_\alpha = \psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} - d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}$,
\begin{align*}
\theta_\alpha \wedge \omega_\alpha &= \big(\psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} - d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}\big) \wedge \psi_{\alpha\beta}\, \omega_\beta \\
&= \psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta - d\psi_{\alpha\beta} \wedge \omega_\beta,
\end{align*}
using $\psi_{\alpha\beta}^{-1} \psi_{\alpha\beta} = I$ in both factors (note that $\psi_{\alpha\beta}$ is a 0-form, so it commutes through the wedge). Summing,
\begin{align*}
(d_A \omega)_\alpha = d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha = d\psi_{\alpha\beta} \wedge \omega_\beta + \psi_{\alpha\beta}\, d\omega_\beta + \psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta - d\psi_{\alpha\beta} \wedge \omega_\beta.
\end{align*}
The $\pm d\psi_{\alpha\beta} \wedge \omega_\beta$ terms cancel, leaving
\begin{align*}
(d_A \omega)_\alpha = \psi_{\alpha\beta}\, d\omega_\beta + \psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta = \psi_{\alpha\beta}\, (d\omega_\beta + \theta_\beta \wedge \omega_\beta) = \psi_{\alpha\beta}\, (d_A \omega)_\beta.
\end{align*}
[guided]
This is the load-bearing computation. We verify that the locally defined candidates $(d_A \omega)_\alpha$ satisfy the same transformation law $\omega \mapsto \psi_{\alpha\beta}\, \omega$ as the local representatives of an $E$-valued form, so they piece together into a global $E$-valued $(i+1)$-form.
Step 1: differentiate the overlap condition $\omega_\alpha = \psi_{\alpha\beta} \omega_\beta$ using the matrix Leibniz rule
\begin{align*}
d\omega_\alpha = d(\psi_{\alpha\beta} \omega_\beta) = d\psi_{\alpha\beta} \wedge \omega_\beta + \psi_{\alpha\beta}\, d\omega_\beta.
\end{align*}
Observe the sign: $\psi_{\alpha\beta}$ is a 0-form so the graded Leibniz sign is $(-1)^0 = +1$. This produces an unwanted cross term $d\psi_{\alpha\beta} \wedge \omega_\beta$ that we must cancel.
Step 2: expand the second contribution using the connection transformation law. Using $\theta_\alpha = \psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} - d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}$ (from Step 1),
\begin{align*}
\theta_\alpha \wedge \omega_\alpha &= \big(\psi_{\alpha\beta}\, \theta_\beta\, \psi_{\alpha\beta}^{-1} - d\psi_{\alpha\beta}\, \psi_{\alpha\beta}^{-1}\big) \wedge (\psi_{\alpha\beta}\, \omega_\beta) \\
&= \psi_{\alpha\beta}\, \theta_\beta\, (\psi_{\alpha\beta}^{-1} \psi_{\alpha\beta}) \wedge \omega_\beta - d\psi_{\alpha\beta}\, (\psi_{\alpha\beta}^{-1}\psi_{\alpha\beta}) \wedge \omega_\beta \\
&= \psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta - d\psi_{\alpha\beta} \wedge \omega_\beta.
\end{align*}
The connection transformation law was designed to produce exactly the term $-d\psi_{\alpha\beta} \wedge \omega_\beta$ needed to cancel the unwanted cross term from Step 1.
Step 3: add the results:
\begin{align*}
(d_A \omega)_\alpha &= d\omega_\alpha + \theta_\alpha \wedge \omega_\alpha \\
&= \big(d\psi_{\alpha\beta} \wedge \omega_\beta + \psi_{\alpha\beta}\, d\omega_\beta\big) + \big(\psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta - d\psi_{\alpha\beta} \wedge \omega_\beta\big) \\
&= \psi_{\alpha\beta}\, d\omega_\beta + \psi_{\alpha\beta}\, \theta_\beta \wedge \omega_\beta \\
&= \psi_{\alpha\beta}\, \big(d\omega_\beta + \theta_\beta \wedge \omega_\beta\big) = \psi_{\alpha\beta}\, (d_A \omega)_\beta.
\end{align*}
The $d\psi_{\alpha\beta} \wedge \omega_\beta$ terms cancel exactly — this is the whole point of the inhomogeneous transformation law for $\theta_\alpha$. If $\theta_\alpha$ transformed as a tensor (i.e. without the $d\psi$-correction), $d + \theta \wedge$ would not define a global operator; the inhomogeneity is precisely what compensates for the failure of $d$ to commute with $\psi_{\alpha\beta}$.
[/guided]
[/step]
[step:Glue the local candidates into a global $E$-valued form and verify linearity]
The compatibility $(d_A \omega)_\alpha = \psi_{\alpha\beta}\, (d_A \omega)_\beta$ on all overlaps is exactly the condition for the local sections $(d_A \omega)_\alpha \in \Gamma(E|_{U_\alpha} \otimes \Lambda^{i+1} T^*U_\alpha)$ to glue into a global section $d_A \omega \in \Omega^{i+1}(E)$ by the sheaf property of sections of a vector bundle. The resulting map
\begin{align*}
d_A: \Omega^i(E) &\to \Omega^{i+1}(E) \\
\omega &\mapsto d_A \omega
\end{align*}
is $\mathbb{R}$-linear because on each chart both $d$ and $\theta_\alpha \wedge (-)$ are $\mathbb{R}$-linear. For $i = 0$, a section $s \in \Gamma(E) = \Omega^0(E)$ has local representative $s_\alpha$ and $(d_A s)_\alpha = ds_\alpha + \theta_\alpha s_\alpha$, recovering the original connection. Naturality — independence of choice of trivialising atlas — follows because the formula is determined by the same data $\{\theta_\alpha\}$ that defines $A$, and the overlap computation shows the output is independent of which chart we use at a point.
[/step]