What you are looking for nowadays goes by the name of the *Rumin complex* and is defined on any contact manifold.  Moreover, there is a vast generalization of this that sometimes goes by the name of 'the variational bicomplex' and sometimes by the name 'characteristic cohomology'.  Here is a brief description that is suited for the question you asked:

On $M = \mathbb{R}^3$ with coordinates $x,y_0,y_1$, consider the differential ideal $\mathcal{I}\subset\Omega^\ast(M)$ generated by the $1$-form $\omega = dy_0 - y_1\ dx$, i.e., $\mathcal{I}$ is the set of linear combinations of all multiples of $\omega$ and $d\omega$.  Note that $\mathcal{I}$ is a homogeneous ideal and equals $\Omega^\ast(M)$ in degrees 2 and 3.  Because $\mathcal{I}$ is closed under exterior derivative, it is a sub-complex of $\bigl(\Omega^\ast(M),d\bigr)$. Thus, there is a graded quotient complex, call it $\bigl(\mathcal{Q},\bar d\bigr)$, that vanishes in degrees above $1$.  Note that $\mathcal{Q}^0 = \Omega^0(M)= C^\infty(M)$, since $\mathcal{I}$ vanishes in degree $0$.

Now, say that an element $\phi \in \mathcal{Q}^1$ is *exact* if $\phi = \bar d f$ for some $f\in \mathcal{Q}^0 = \Omega^0(M)= C^\infty(M)$.  Unfortunately, unlike $\bigl(\Omega^\ast(M),d\bigr)$, the complex $\bigl(\mathcal{Q},\bar d\bigr)$ is not locally exact in positive degree.  In fact, $\bar d \phi =0$ for all $\phi\in\mathcal{Q}^1$, even though $\bar d: \mathcal{Q}^0\to \mathcal{Q}^1$ is not onto.

Let me pause just a moment to explain how this fits into your question.  Your equation $P(x)y'' + Q(x)y' +R(x)y = 0$ is encoded as the $1$-form $\phi = P(x) dy_1 + (Q(x)y_1 + R(x) y_0) dx$ (which represents the same class as the $1$-form $P(x) dy_1 + Q(x)dy_0 + R(x) y_0 dx$ in $\mathcal{Q}^1$), and you are asking when there is a function $f(x,y_0,y_1)$ such that $\phi = \bar d f$.  (You should verify that $f = P(x) y_1 + (Q(x)-P'(x))y_0$ works when your equation is satisfied and that, otherwise, nothing does.)

Now, how can we test for exactness in this sense?  This is where the Rumin complex (aka the variational bicomplex, etc.) comes in.  It turns out that there is a way to embed the operator $\bar d:\mathcal{Q}^0\to\mathcal{Q}^1$ into a complex that provides a fine resolution of the constant sheaf, the same way that the exterior derivative does for the full complex of exterior differential forms.  

What you do is this:  Let $\mathcal{E}^2\subset\Omega^2(M)$ be the set of multiples of $\omega$ and let $\mathcal{E}^3=\Omega^3(M)$.  

We now want to define a complex
$$
0\longrightarrow \mathcal{Q}^0
\buildrel{\bar d}\over\longrightarrow \mathcal{Q}^1
\buildrel{D}\over\longrightarrow \mathcal{E}^2
\buildrel{d}\over\longrightarrow \mathcal{E}^3
\longrightarrow 0.
$$
The map from $\mathcal{E}^2$ to $\mathcal{E}^3$ is the usual exterior derivative, so the only thing left to define is the map $D:\mathcal{Q}^1\to \mathcal{E}^2$.  To do this, we first define a (first-order) operator $\delta: \mathcal{Q}^1\to\Omega^1(M)$, by requiring that $\delta(\phi)$ be a 1-form representing $\phi$ in the quotient complex and that $d\bigl(\delta\phi\bigr)$ lie in $\mathcal{E}^2$, i.e., that it be a multiple of $\omega$.  (I'll let you write down the formula for $\delta$ in local coordinates.)  Now, define $D\phi$ to be $d\bigl(\delta\phi\bigr)$.  (Sounds almost trivial doesn't it?)  The operator $D$ is easily verified to be second order and linear.

Now, it is not hard to verity that this complex is locally exact in positive degrees.  (It also gives a fine resolution of the constant sheaf, so its cohomology on $M$ is canonically isomorphic to the deRham cohomology of $M$.)  In particular, the local condition that $\phi\in\mathcal{Q}^1$ be exact is that $D\phi=0$.  

You should verify (after you have defined $D$) that this reproduces your condition precisely in the linear case you asked about.