Every smooth function on a manifold defines a hamiltonian vector field. Let $(M,\omega)$ be a symplectic manifold, and let $H$ be a smooth function on $M$. I want to show that $H$ is a Hamiltonian function i.e. there exists a smooth vector field $X$ on $M$ such that
$$\iota_X\omega=dH \text{ }(*)$$
Since $\omega$ is a non-degenerate form as it's a symplectic form, then we can see that we can solve $(*)$ for $X$.
Question: Why can we do this?
My idea was to expand this equation in local coordinates. For example, let's do a simple example when $\dim(M)=2$. Let's choose some point $p\in M$ and some chart $(U,\varphi)$ containing that point with $\varphi(p)=0$. Also, consider some symplectic form $\omega_p=c(x,y)(dx)_p\wedge (dy)_p$ where $c(x,y)$ is non-zero (I am going to drop a subscript $p$ for a simplification).
For a given function $H:M\to\mathbb{R}$, we want to find $X=a(x,y)\partial_x+b(x,y)\partial_y$ where $\partial_x=\frac{\partial}{\partial x}$ s.t. $(*)$ holds in local coordinates i.e. we want to solve it for the functions $a(x,y)$ and $b(x,y)$.
Since $\iota_X\omega$ and $dH$ are linear, then it's enough to check $(*)$ on the basis $\{\partial_x,\partial_y\}$ of $T_pM$. I will write $H_x$ instead of $\frac{\partial H}{\partial x}$. Then we can see that $\iota_X\omega(\partial_x)=dH(\partial_x)$ gives us
$$dH(\partial_x)=(H_xdx+H_ydy)(\partial_x)=H_x\text{ and }$$
$$\iota_X\omega(\partial_x)=\omega(X,\partial_x)=c(x,y)dx\wedge dy(a(x,y)\partial_x+b(x,y)\partial_y,\partial_x)=-c(x,y)b(x,y)$$
So, we have that
$$H_x=-c(x,y)b(x,y)$$
The same argument works for $\partial_y$ and gives us
$$H_y=c(x,y)a(x,y)$$
Since $c(x,y)\neq0$ then we have that
$$a(x,y)=\frac{H_y}{c(x,y)}\text{ and }b(x,y)=\frac{-H_x}{c(x,y)}$$
which are smooth function. Therefore, $X$ is a smooth vector field which satisfies $(*)$.
So, as I understand, I can use the same approach for the bigger dimension. Where instead of $c(x,y)\neq0$, I will use the fact that $\omega$ is non-degenerate? So, in other words, if I have $H$ and $\omega$, then I can explicitly find coordinates of $X$ by solving a similar system.
 A: The proof that you can "solve for $X$ "is actually a rather simple extension of a linear algebra fact, which is why I'll treat the vector space case closely. Let's recall what non-degeneracy means.

Definition.
Let $V$ be a finite-dimensional vector space over $\Bbb{R}$, and let $\omega:V \times V \to \Bbb{R}$ be bilinear (and skew-symmetric... but this isn't really necessary). We say $\omega$ is non-degenerate if the map $\omega^{\flat}:V \to V^*$ defined by
\begin{align}
\omega^{\flat}(x):= \omega(x, \cdot) \equiv \iota_x \omega \equiv \bigg( y \in V \mapsto \omega(x,y) \in \Bbb{R}\bigg) \in V^*
\end{align}
is injective (or in finite-dimensions, we can equivalently require that it be an isomorphism). Note that I use $\equiv$ to mean "same thing expressed in different notation".

You may have seen the definition probably stated as "for all $x \in V$, if for all $y \in V$, $\omega(x,y) = 0$ then $x=0$". Well, this is exactly what it means for $\omega^{\flat}$ to be injective (and hence an isomorphism).
Now, being an isomorphism means it has a linear inverse, which we may denote as $\omega^{\sharp}:V^* \to V$. So, for any covector $\alpha \in V^*$, we can consider the vector $x:= \omega^{\sharp}(\alpha) \in V$. What's special about this vector $x$? Well, just apply $\omega^{\flat}$ to both sides of this equation and you'll see that
\begin{align}
\omega^{\flat}(x) = \omega^{\flat}(\omega^{\sharp}(\alpha)) = \alpha
\end{align}
in other words,
\begin{align}
\omega(x, \cdot) = \iota_x\omega = \alpha
\end{align}
This is why given a covector $\alpha$, we can always find a vector to make the above equation true.

In your case, you just have to repeat everything pointwise. $dH$ is a covector-field (i.e a $1$-form). So, consider the vector field $X$ defined pointwise as $X_p := (\omega_p)^{\sharp}\left( dH_p\right) \in T_pM$. Then, it will satisfy
\begin{align}
\omega_p(X_p, \cdot) = dH_p
\end{align}
i.e if you remove the point $p$, then $\iota_X \omega = \omega(X, \cdot) = dH$.

By the way, doing thing in coordinates may be a little hard, because as you can see, it involves the inverse mapping $\omega^{\sharp}$. But anyway, if you're working in some chart $(U,x)$ of the manifold $M$, with the coordinate basis $\{\partial/\partial x^1, \dots \partial/ \partial x^n\}$, and dual basis $\{dx^1, \dots dx^n\}$, then define the functions
\begin{align}
\omega_{ij}:= \omega\left( \dfrac{\partial}{\partial x^i}, \dfrac{\partial}{\partial x^j}\right)
\end{align}
and let $[\omega^{ij}]$ be the inverse matrix of $[\omega_{ij}]$. Then, the components of the vector field $X = \sum_{i}X^i\frac{\partial}{\partial x^i} $ will be
\begin{align}
X^i &= \sum_{i=1}^n\omega^{ij} \dfrac{\partial H}{\partial x^j}.
\end{align}
(you see, the appearance of the inverse matrix entries makes things not-so-easy)
