I'm trying to understand the proof of of Proposition 2.5.3 in $\textit{An Introduction to Mechanics and Symmetry}$ by Marsden and Ratiu. They claim that a vector field $X$ on a symplectic, simply connected manifold $(M, \omega)$ (say ($T^*(\mathbb{R^n}), \omega)$) is Hamiltonian if and only if the Jacobian matrix of $X$ is $\omega$-skew, that is, if $DX \in \frak{sp}_n(\mathbb{R})$ (in canonical coordinates).
If we set $\theta = \iota_X(\omega)$, then $\theta$ is closed if and only if $d\iota_X(\omega) = 0$.
Now, $\mathcal{L}_X(\omega) = \iota_Xd\omega + d\iota_X(\omega)$ by Cartan's Magic Formula. Since $\omega$ is closed, $d\omega = 0$, so $\iota_Xd\omega = 0$.
Now, by another formula, $\mathcal{L}_X(\omega)(Y,Z) = X(\omega(Y,Z)) - \omega(\mathcal{L}_X(Y),Z) - \omega(Y,\mathcal{L}_X(Z))$.
I'm assuming $X(\omega(YZ)) = \iota_Xd\omega$.
In that case, $-d\iota_X(\omega) = \omega(\mathcal{L}_X(Y),Z) + \omega(Y,\mathcal{L}_X(Z))$.
I'm further assuming $\omega(\mathcal{L}_X(Y),Z) + \omega(Y,\mathcal{L}_X(Z)) = \omega(DX \cdot Y,Z) + \omega(Y, DX \cdot Z)$.
If both of these assumptions are true, then $\theta$ is closed if and only if $DX \in \frak{sp}_n(\mathbb{R})$, because $\omega(DX \cdot Y,Z) = -\omega(Y, DX \cdot Z)$ then, in which case $\theta$ is exact since $M$ is simply connected, so $\theta = \iota_X(\omega) = dH$, as desired.
Can anyone help me out with verifying/correcting the assumptions?
If you're in local canonical coordinates, then the proof is exactly the same as the linear case, and is given in Proposition 2.5.1. I understand you're trying to come up with a "coordinate-free" proof, but the problem is that a quantity like $DX$ has no invariant meaning, and is only defined with respect to a coordinate system.
For a proof more along the lines you're attempting, imagine we have a canonical coordinate system $(z^1, \ldots, z^{2n})$ (so $\left[\omega\left(\frac{\partial}{\partial z^i},\frac{\partial}{\partial z^j}\right)\right] = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n\end{bmatrix}$). As you rightly say, $\mathcal{L}_X\omega = d(i_X\omega)$. However, if $X = X_H$ is a Hamiltonian vector field, then $i_{X_H}\omega = dH$, and so $\mathcal{L}_X\omega = d(dH) = 0$. If follows that $$ X(\omega(Y,Z)) = \omega(\mathcal{L}_XY,Z) + \omega(Y,\mathcal{L}_X Z). $$ Now take $Y = \frac{\partial}{\partial z^i}$ and $Z=\frac{\partial}{\partial z^j}$. The left hand side will vanish, since the components of $\omega$ are constant in the symplectic basis. Also $$ \mathcal{L}_XY = \left[X,\frac{\partial}{\partial z^i}\right] = -\frac{\partial X^m}{\partial z^i}\frac{\partial}{\partial z^m} = -(DX)_i^{\ m}\frac{\partial}{\partial z^m}, $$ and similarly for $\mathcal{L}_XZ$. This will give the result you want.