Suppose that I have a small family of partial differential equations such that $y(x)$ must be a solution to $$ P_i(x,D) y(x)=0,\ \ \ \ \ \ i=1,...,m $$ for all $i=1,...,m$, where I am using multi-index notation $x=(x_1,...,x_n)$, and as well $P_i(x,D)=\sum_{j=0} p_{ij}(x)(\frac{\partial}{\partial x_1})^{\alpha_1}\cdots (\frac{\partial}{\partial x_n})^{\alpha_n}$, and $p_{ij}(x)\in \mathbb C[x_1,...,x_n]$. My question is whether there exists a fast algorithm for determining if $P_i$ describe a holonomic D-module.
I am reading about Hilbert polynomials, rings, modules, and characteristic varieties but it would be nice if there were an elementary method, at least for a special system, e.g. if n=m and all operators $P_i$ depend only on $x_i$ and $\frac{\partial}{\partial x_i}$. Is this at all true? Or must I calculate the Hilbert polynomial everytime?
Let $X$ be a scheme locally of finite type over $\mathbf{C}$. You can analytify it : you associate to $X$ a complex analytic space $X^{\textrm{an}}$. If $X$ is a closed subscheme of an affine space, defined by polynomials $Q_1,\ldots,Q_d$, the analytic space $X^{\textrm{an}}$ will be the closed analytic subspace of the affine space defined by the same polynomials. In fact, the underlying topological space of $X^{\textrm{an}}$ is equal to the set $X(\mathbf{C})$ of closed points of $X$ in this case, and you can glue this to get the general construction. This gives you a functor $X\mapsto X^{\textrm{an}}$, which has very nice properties. This is GAGA theorem. (You can see J.-P. Serre's seminal paper Géométrie algébrique et géométrie analytique on the subject, but also these Zhao's notes for a modern treatment.) This functor concerns certains modules over $\mathscr{O}_X$-modules, and tells that their "image" by this functor has the same properties than the properties that they have. Now, $\mathscr{D}_{X/\mathbf{C}}$ makes sense (see EGA IV.4 paragraph 16, section 8, definition (16.8.1)) and if $X$ is smooth over $\mathbf{C}$, holonomicity makes sense and has all expected nice properties. Now, one can also show that GAGA can be extended to the case of $\mathscr{D}_{X/\mathbf{C}}$-modules (left, let's say) and that if $\mathscr{M}$ is a (left) $\mathscr{D}_{X/\mathbf{C}}$-module of dimension $d$, so is its analytification $\mathscr{M}^{\textrm{an}}$. Thus our functor sends holonomics to holonomics.
Why all of this ? Just to say that even if you are looking to analytic solutions to your complex differential system, as it is defined by polynomial equations, you can calculate its dimension and test its holonomicity in the algebraic case. You are in the case $X = \mathbf{A}_{\mathbf{C}}^n$ of the affine space of dimension $n$ over $\mathbf{C}$ so that you can work with global sections, and for this, yes, you have an algorithm.
You have Macaulay 2, which is pretty well-known in the effective world of algebraic geometry and commutative algebra, and it has a module dedicated to $\mathscr{D}$-modules, with a documentation in post-script here, that redirects to a bibliography of papers from which algorithms have been used. I used this a long time ago to calculate $b$-functions of $\mathscr{D}$-modules, and it was working pretty well. In your special case you will want to calculate the dimension of your $\mathscr{D}$-module and test if it is holonomic by checking it this dimension is equal to $n$ or not. The function Ddim does this job.