A complex function $f$ (domain and image a subset of $\mathbb{C}$; I'm intentionally vague here because this is part of my confusion) is called holonomic if it satisfies a differential equation with polynomial coefficients, i.e. there are $q_0(z),\cdots,q_d(z)\in\mathbb{C}[z]$ not all zero such that $$ q_0(z)f(z)+q_1(z)f'(z)+\cdots+q_d(z)f^{(d)}(z)=0 $$ It is claimed often in literature that "holonomic functions can only have finitely many singularities". I'm curious as to the precise meaning of this statement, and a proof of the statement itself.
For instance, on the wikipedia page about holonomic functions, it is claimed that the function $g$ given by $z\mapsto \frac{z}{e^z-1}$ is not holonomic because it has infinitely many singularities. Now, in this particular case $g$ is meromorphic, and I can see how the statement makes sense (but still don't know how to prove it for general meromorphic functions).
For instance I can see how a holonomic meromorphic function $f$ has finitely many poles, because if it has a pole at $z_0$ of order $m$ then (I think) $f^{(j)}$ has a pole of order $m+j$ at $z_0$, hence by $$ f^{(d)}(z)=\frac{q_0(z)f(z)+q_1(z)f'(z)+\cdots+q_{d-1}(z)f^{(d-1)}(z)}{q_d(z)} $$ we find that $z_0$ must be a zero of $q_d$. However this still doesn't cover the case of essential singularities.
It gets even trickier in the following situation: Suppose that $h$ is given by a power series that is convergent within the open unit disk $D$, but not at a dense subset of its boundary. Thus $D$ forms a natural boundary for $h$ (see http://mathworld.wolfram.com/NaturalBoundary.html). Can we conclude from the above "statement" that $h$ is not a holonomic function? What about the power series of $h$?