Codimension of the range of certain linear operators Assume that $P(x,y), Q(x,y) \in \mathbb{R}[x,y]$  are two polynomials. We define a linear map
$D$ on $\mathbb{R}[x,y]$ with $D(U)=PU_{x}+QU_{y}$. In fact $D$ is the differential operator corresponding to the vector field $P\partial_{x}+Q\partial_{y}$.

Are there  polynomials $P$ and $Q$ such that the codimension of the range of $D$ is finite but different from $0$ and $1$?

Note  This codimension is $0$, $1$ and $\infty$ for $\partial_{x}$,    $x\partial_{x}+y\partial_{y}$  and  $x\partial_{x}-y\partial_{y}$, respectively.
Motivation It can be easily shown that this codimension is an upper bound for the number of closed orbits of the polynomial vector field $P\partial_{x}+Q\partial_{y}$.
However there is  a  Painful fact: If  an  smooth  vector  field  $X$ on $\mathbb{R}^{2}$ has  a  limit  cycle which  surrounds a  nondegenerate  singularity (or  a  singularity which  has  a  local smooth first integral at  a  deleted  neighborhood of the  singularity and is  discontinuous at the  singularity), then  the  codimension of the  range  of the  differential  operator  $D_{X}$ on $C^{\infty} (\mathbb{R}^{2})$  is  infinite ! The reason is as follows.
Assume  that  a  limit  cycle  $\gamma$  is  attractor  and surounds  a  source singularity  at the  origin.
Non degeneracy of the singularity help  us  to  find an open set  D  around  origin whose  boundary  $S$ is a smooth  closed  curve and the vector field is transverse to $S$ toward the  exterior  of  $S$.  So the  flow of  $X$ defines  a smooth  retraction  $r:\bar{D}\setminus \{0\} \to S$. In fact $r(x)$ is the intersection point  of  the  orbit of $x$  with $S$  This  shows  that there  is  an  smooth real valued  function $\phi$  on  $\mathbb{R}^2 \setminus  \{0\}$  such that locally  around the  origin we  have $X. \phi =0$   and $P\circ \phi$  is  discontinuous at origin for  every  polynomial $P(x)$. Such  $\phi$ can  be  constructed locally  by  $\psi \circ 
r$ where $\psi:S \to \mathbb{R}$ is an arbitrary non constant smooth function.Now we choose  an arbitrary smooth extension of this locally defined  $\phi$ to whole punctured plane. This extension is denoted again
by $\phi$. Obviously every nontrivial linear combination $\sum \lambda_i \phi^i$ is discontinuos at origin.  Because its  restriction to  every  small  neighborhood  of the  origin has  the  same  range(image) as the  range  of  $\sum \lambda_i \psi^i: S \to  \mathbb{R}$. However $X.\phi^i$ is a smooth function on whole plane, for every $i$  since  it  vanish  locally  around the  origin.
Now for every $n$,the following set represents an independent subset of  $C^{\infty}(\mathbb{R}^2)/Range(D_X)$
$$\{X.\phi, X.\phi^2, \ldots, X.\phi^n \}$$  To prove this, we  strongly  use  the  fact that  we  have  at  least  one  limit  cycle.
More precisely: assume that $\sum_{i=1}^n  \lambda_i  (X.\phi^i)= X.F$ for some smooth function $F\in C^{\infty} (\mathbb{R}^2)$. Then $F-\sum \lambda_i \phi^i$ is a first integral(constant of motion) on the puntured plane. So it is constant on the basian of attraction of the limit cycle. A punctured neighborhood of the origin is contained in the basian of attraction of limit cycle. Thus locally around the origin,$F$ differs $\sum_{i=1}^n  \lambda_i \phi^i$ by a constant(In a deleted   neighborhood). This is a   contradiction because  $F$ is  smooth  at origin  but $\sum  \lambda_i \phi^i$ is  discontinuous at origin.
So  unfortunately we encounter with the following  two inconsistent and opposite but true results:

1)The  codimension  of  the  range  of  $D_X$ is an upper bound for the number  of limit  cycles.


2.If the codimension is finite then there is no any limit cycle.

So  we  should  choose  an  appropriate  function  algebra, different from $C^{\infty} (\mathbb{R}^{2})$,  which is  invariant under the  differential operator corresponding  to  an  algebraic  vector  field.
#Added:
According to the answer of  Loic  Teyssier we put $X=x^k\partial_x+y\partial_y$.
Then $X$ can act on various  spaces $C^{\infty}(\mathbb{R}^2),\; C^{\omega}(\mathbb{R}^2)$ or the space  of  holomorphic functions  from $\mathbb{C}^2$ to $\mathbb{C}$. What can be said  about the  dimension of  the  cokernels in each of these actions?
 A: The question admits a positive answer when one looks at it for the ring of formal power series $\mathbb R[[x,y]]$. For instance the vector field $x^k\partial_x+y\partial_y$ has a formal cokernel of dimension $k$. Yet when restricted to polynomials the property disappears: the polynomial cokernel is infinite-dimensional (you can never reach monomials of the form $x^{k-1}y^m$ for $m\in\mathbb N$ nor $x^n$ with $n<k$, and those are the only problems).
Somehow I believe the answer should also be positive in the polynomial case. But I have no example offhand.
A: This is not an answer. 
I've tried some experiments in Maple and I believe that the answer to your question is negative. I've looked at the infinite matrix given by $D$ as I've suggested in my comments.
A column in this infinite matrix is given by coefficients of $D(x^iy^j) = iP(x,y)x^{i-1}y^j + jQ(x,y)x^iy^{j-1}$ which means that the entries in the matrix look like $\lambda p_{ab} + \mu q_{cd}$ for integral $\lambda$ and $\mu$. If one such an entry is zero, then it implies vanishing of infinitely many entries in the matrix, which are given by multiples of this specific $(\lambda, \mu)$ coming from some higher powers $x^{ki}y^{kj}$.
For the corner cases $x^i$ and $y^j$ one recovers multiples $P$ and $Q$, the coefficients are just more spread out throughout the rows. Since the matrix is band diagonal (the girth given by degrees of $P$ and $Q$), it seems that there is a slim that the codimension could be finite and bigger than 1.

edit:
If you want to consider various function spaces then it may be fruitful to reinterpret the statement in the theory of $\mathcal{D}$-modules, see for  example chapter $8$ of "Homological Algebra"  by S. I.  gelfand  and Yu. I  Manian. It will also give you more literature to study the broader scope of your problem. I've took the liberty to add appropriate tags to your question. Please take the following lines with a grain of salt as it has been quite long since I've dealt with $\mathcal{D}$-modules.
Let $\mathcal{O}$ denote the space of polynomials in two variables and let $\mathcal{D}$ be the corresponding Weyl algebra, i.e. the space of differential operators in two variables with coefficients from $\mathcal{O}$. Finally, let $M =\mathcal{D}/\mathcal{D}D$ denote the quotient of $\mathcal{D}$ by the left ideal generated by $D$. If you apply the left exact hom-functor $\mathrm{Hom}(_,\mathcal{O})$ on the short exact sequence
$$
0 \to \mathcal{D} \xrightarrow{Q \mapsto QD} \mathcal{D} \to M \to 0
$$
you will obtain the following long exact sequence
$$
0 \to \mathrm{Hom}_\mathcal{D}(M,\mathcal{O}) \to \mathrm{Hom}_\mathcal{D}(\mathcal{D},\mathcal{O}) \to \mathrm{Hom}_\mathcal{D}(\mathcal{D},\mathcal{O}) \to \mathrm{Ext}^1_\mathcal{D}(M,\mathcal{O}) \to \mathrm{Ext}^1_\mathcal{D}(\mathcal{D},\mathcal{O}) \to \ldots
$$
The vector space $\mathrm{Hom}_\mathcal{D}(\mathcal{D},\mathcal{O})$ is actually isomorphic to $\mathcal{O}$ via $\varphi \mapsto \varphi(1)$ and if you unwind the definition of the hom-functor and use this isomorphism, you'll see that the first space $\mathrm{Hom}_\mathcal{D}(M,\mathcal{O})$ is actually isomorphic to $\mathrm{Ker}(D)$ -- the kernel of $D$ acting on the space of polynomials $\mathcal{O}$. Since $\mathrm{Ext}^1_\mathcal{D}(\mathcal{D},\mathcal{O}) = 0$ we obtain exact sequence
$$
0 \to \mathrm{Ker}(D) \to \mathcal{O} \xrightarrow{f \mapsto D(f)} \mathcal{O} \to \mathrm{Ext}^1_\mathcal{D}(M,\mathcal{O}) \to 0,
$$
which tells us that $\mathrm{Ext}^1_\mathcal{D}(M,\mathcal{O})$ is actually the cokernel of $D$. All of this goes through for a general $\mathcal{D}$-module $N$ instead of $\mathcal{O}$, the only change is that the kernel is now taken in the space $N$. Thus your question can be reinterpreted/generalized as:
Is there a $\mathcal{D}$-module $N$ such that $\mathrm{Ext}^1_\mathcal{D}(\mathcal{D}/\mathcal{D}D,N)$ is finite-dimensional of dimension greater that one?
There are several computer packages that can give you the cokernel (or ext) of any differential operator over a Weyl algebra. I wouldn't bet on them being able to handle any "nonalgebraic" modules like the space of real analytic functions, but you can try to experiment with some extensions of $\mathcal{O}$.
There are two articles by Coutinho that investigate on the topic: Extensions of modules over Weyl algebras and Modules of codimension one over Weyl algebras).
It would seem that the corollary 2.6 of the former article answers your question in negative since $\mathcal{O} = \mathcal{D}/(\mathcal{D}\partial_x + \mathcal{D}\partial_y)$ is a holonomic module and thus it is singular. But I have no idea whether there can exists $D$ such that $\mathcal{D}/\mathcal{D}D$ is nonsingular.
