Coordinate free proof that curvature is the "square" of the connection

Here's the setup. Consider a vector bundle $E$ over a manifold $M$ and let $\Omega^*(M, E)$ denote the space of $E$-valued differential forms (i.e. the space of sections of the vector bundle $\bigwedge^* T^*M \otimes E$). Let $\nabla$ be an affine connection on $E$, i.e. a map $\Omega^0(M, E) \to \Omega^1(M, E)$ satisfying the Leibniz rule. For each vector field $X$ there is a contraction operator $\iota(X)$ on $\Omega^*(M,E)$, and we form $\nabla_X = \iota(X) \nabla$.

Recall that the curvature of $\nabla$ is defined to be $F(X,Y) = [\nabla_X, \nabla_Y] - \nabla_{[X,Y]}$. It is well known that $F$ is actually a $End(E)$-valued 2-form. $\Omega^2(M, End(E))$ acts on $\Omega^*(M, E)$ in the standard way, and it is also well known that $\nabla^2: \Omega^k(M,E) \to \Omega^{k+2}(M, E)$ is given by the action of $F$.

My question is about the proof that "$\nabla^2 = F$" in the sense described above. The textbook proofs all involve trivializing $E$ and doing a coordinate calculation; I don't mind trivializing $E$, but for my own nefarious (and perhaps misguided) purposes I would like to avoid coordinates. Here is what I mean.

If $E$ is trivial then $\nabla = d + \omega$ where $\omega$ is an $End(E)$-valued 1-form. So $F(X,Y) = [\iota(X)(d+\omega), \iota(Y)(d + \omega] - \iota([X,Y])(d + \omega)$, and it seems to me that one should be able to carry out the argument using only Cartan's homotopy formulas. There are complications which arise from the fact that the $d$ operator and the contraction operators are different for $\Omega^*(M, End(E))$ versus $\Omega^*(M,E)$, but nevertheless I can get tantalizingly close. I can show my progress if it would help you help me.

In any event, I've been at it for quite some time now and I've decided it's time to seek assistance. So I'm wondering if anybody knows either how to do this, a reason why it can't be done, or a place where I could look it up. Thanks in advance!

Let $\sigma \in \Omega^n(M,E)$ and $X$, $Y$, $Z_1$, . . ., $Z_n$ be vector fields, then:

$(\nabla^2\sigma)(X,Y,Z_1,. . . , Z_n) = \nabla_X(\nabla(\sigma))(Y,Z_1,. . . , Z_n) $

$= \nabla_X(\nabla_Y(\sigma))(Z_1,. . . , Z_n) -\nabla\sigma(d_X(Y))(Z_1,. . . , Z_n) $

$= \nabla_X(\nabla_Y(\sigma))(Z_1,. . . , Z_n) -(\nabla_{d_X(Y)}\sigma)(Z_1,. . . , Z_n) $

$= -\nabla_Y(\nabla_X(\sigma))(Z_1,. . . , Z_n) +(\nabla_{d_Y(X)}\sigma)(Z_1,. . . , Z_n) $

$= \frac{1}{2}(\nabla_X(\nabla_Y(\sigma)) -\nabla_Y(\nabla_X(\sigma)))(Z_1,. . . , Z_n) -\frac{1}{2}(\nabla_{d_X(Y)-d_Y(X)}\sigma)(Z_1,. . . , Z_n) $

$= \frac{1}{2}(([\nabla_X, \nabla_Y])(\sigma))(Z_1,. . . , Z_n) -\frac{1}{2}(\nabla_{[X,Y]}\sigma)(Z_1,. . . , Z_n) $