I will start commenting on Mariano's answer. I believe it is a perfect answer for the question

How do sheaves arise in studying
solutions of differential equations ?

but not for the question

How do sheaves arise in studying
solutions to **ordinary** differential
equations ?

According to the current terminology a function $f$ satisfying $X(f)=0$ is not a solution of the vector field $X$ but a first integral. Moreover, if $X = a(x,y) \partial_x + b(x,y) \partial_y$ then
$$
X(f) = a \partial_x f + b \partial_y f .
$$
Thus $X(f)=0$ is a PDE and not an ODE. Indeed t3suji made the same point at a comment on Mariano's answer. I understand the solutions of (the ODE determined by) $X$ as functions $\gamma : V \subset \mathbb R \to U$ satisfying $X(\gamma(t))=\gamma'(t)$ for every $t \in V$. Notice that here indeed we have a system of ODEs.

A vector field can be thought as autonomous differential equation and I do not see clearly how to consider the sheaf of its solutions.

On the other hand when we have a non-autonomous ordinary differential equation then there is its sheaf of solutions. This sheaf is a sheaf over the time variable
only and not the whole space. ( At this point it is natural to talk about connections and/or jet bundles but I will try to keep things as elementary as possible. )

Note that in general the sheaf of solutions will not be a sheaf of vector spaces: the sum of two solutions, or the multiplication of a solution by a constant need not to be a solution. This will occur only when the differential equation is linear.

The differential equations $y'(t) = y$ and $y'(t) = y^2$, both defined over the whole real line, are examples of differential equations with non-isomorphic sheaves of solutions. The solutions of the first ODE are the multiples of $\exp t $ and define a sheaf of $\mathbb R$-modules. The solutions of the second ODE are zero and $\frac{1}{\lambda - t}$ with $ \lambda \in \mathbb R$. They do define a sheaf of sets, but not a sheaf of $\mathbb R$-modules.

To obtain examples of linear differential equations with non-isomorphic sheaves, one has to have nontrivial fundamental group on the time-variable of the differential equation. Thus it is natural to consider complex differential equations over $\mathbb C^{\ast}$.

The equations $y'(z) = \frac{ \lambda y(z)}{z}$ parametrized by $\lambda \in \mathbb C$ have non-isomorphic sheaves of solutions. More precisely,

- if $\lambda \in \mathbb Z$ then the solution sheaf is the free $\mathbb C$-sheaf of rank one (solutions of the ODE are complex multiples of $z^{ \lambda }$);
- if $\lambda \in \mathbb Q - \mathbb Z $ then the solution sheaf has no global sections but some tensor power of it does;
- if $\lambda \in \mathbb C - \mathbb Q$ then the solution sheaf has no global sections nor any of its powers does.