Define a sequence of polynomials: $p\_0(x)=1$, $p\_1(x)=2+x$, and, for $n\ge 1$, $p\_{n+1}=(4n+2)p\_n(x)+x^2p\_{n-1}(x)$ so that the first few are $1, x+2, x^2+6x+12, x^2+12x^2+60x+120$. 


Is there an elementary explanation for why $p\_n(x)-e^xp\_n(-x)=O(x^{2n+1})$? 


I understand that the polynomials $\{p\_n\}$ form a variant of the ``reverse Bessel polynomials”, and that the result I seek is tantamount to the fact that $p\_n(x)/p\_n(-x)$ is the (n,n) Pade approximation of $e^x$. I wish to avoid, for expository reasons, all the theories of continued fractions, orthogonal polynomials, Pade approximations, etc. I wouldn’t mind using the fact that $f\_n(x):=e^{-x/2}p\_n(x)$ satisfies the differential equation $4xf\_n’’-8nf\_n’-xf\_n=0$ if that turns out to be helpful. 



Approximation of $e^x$ by rational functions
We can indeed deduce this from the ODE you gave. We want to show that
$$
g(x):=f\_n(x)-f\_n(-x) = O(x^{2n+1})
$$
near $x=0$. First of all, observe that $h(x):=f\_n(-x)$ satisfies the same ODE, and thus so does $g=f\_n-h$. The equation has a regular singular point at $x=0$, and can be discussed with [standard methods:](http://mathworld.wolfram.com/FrobeniusMethod.html) Look for solutions of the form $g=\sum\_{j\ge 0} a\_j x^{j+\lambda}$. The indicial equation is
$$
4\lambda (\lambda-1) - 8n\lambda = 0 ,
$$
which has the solutions $\lambda\_1=2n+1$, $\lambda\_2=0$. So the general solution is of the form $AP(x) + B x^{2n+1}Q(x)$, where $P,Q$ are power series with $P(0),Q(0)\not= 0$. Since $g(0)=0$, the solution we're interested in is of the form $g=Bx^{2n+1}Q(x)$, as claimed.


