How to create a Pade approximation for a difficult function with a divergent Taylor series?

For small values of $a$, the Pade approximants $P^N_M$ calculated from the terms of your power series give a good approximation. Here is a plot of the (log of) ratio of some diagonal Pade approximants/numeric:

This is more of a back-of-envelope calculation: you should play around with it and also look at the $P^N_{N+1}$.

For large values of $a$, split the integral up

$$\tag{1} 2 \sqrt{2\pi} \ h(a)=\int dx \ e^{-(x+a)^2/2}\ln(p) \ +\int dx \ e^{-(x-a)^2/2}\ln(p) $$

Looking at just the first term on the RHS for now, change integration variables $u=x+a$, and expand $\ln (p)$ near the maxima at $u=0$ to get

$$\tag{2} \int du \ e^{-u^2/2}\left[p_0+u^2p_2+\cdots \right] $$

Only terms even in $u$ will survive the integral. The $p_i(a)$ are coefficients of expanding $\ln(p)$. Notice that in the second term (1) we may change variables $u=x-a$ to get an expression identical to (2), so we have

$$\tag{3} h(a)\sim \frac{1}{\sqrt{2\pi}}\int du \ e^{-u^2/2}\left[p_0+u^2p_2+\cdots \right] \qquad, \qquad a \to \infty $$

Term by term integration (thanks Mathematica!) yields

$$\tag{4} h(a)\sim \frac{\ln(8\pi)}{2}-\ln(1+e^{-2a^2})+\frac{1}{2}(1-a^2\operatorname{sech}^2(a^2)) \qquad ,\qquad a\to \infty $$

Here is a plot of the ratio approx/numeric. Note the vertical scale:

It turns out the approximation (4) is quite good even for 'less large' $a$: it has a maximum error of just $6\%$ near $a=3/2$. Of course for small $a$ the power series or Pade approximants are much better.

Using Mathematica:

p[x_, a_] = Map[a |-> PDF[NormalDistribution[a, 1], x], {-a, a}] // Mean // Simplify
h[a_] := NIntegrate[-p[x, a] Log[p[x, a]], {x, -Infinity, Infinity}]
v[x_, n_] := D[-p[x, a] Log[p[x, a]], {a, n}] /. a -> 0
c = Map[a |-> (x^a/a!) Integrate[v[x, a], {x, -Infinity, Infinity}], Range[0, 20, 2]] // Total
f[x_] = PadeApproximant[c, {x, 0, 4}]
Plot[{f[a], h[a]}, {a, -2, 2}]

I get some plot of the function between -2 and +2.

Here is the actual function: $$\frac{(1/2 (1 + \log(2 \pi)) + 1/2 x^2 (4 + 3 \log(2 \pi)) + 1/12 x^4 (22 + 7 \log(2 \pi)))}{(1 + 3 x^2 + (7 x^4)/6)}$$