Step 1: analyze the recurrence
I implemented the recurrence in Python, using basic numpy and scipy.special for the erf
function. Why? Because it is simple and cheap (for reasonable $p$).
import numpy as np
import scipy.special as sp
def F1(z):
return np.sqrt(np.pi/2)*(1.0+sp.erf(z/np.sqrt(2)))
def F2(z):
return F1(z)*z + np.exp(-z**2/2.0)
def F(p,z):
if (p==1):
return F1(z)
elif (p==2):
return F2(z)
else:
fpm2 = F1(z)
fpm1 = F2(z)
for pi in range(3,p+1):
fp = (pi-2)*fpm2 + z* fpm1
fpm2 = fpm1
fpm1 = fp
return fp
I compared results from this to fixed precision results using Maple 2016 (Digits set to 500). The following graph shows the relative error between the recurrence results and Maple results.

You can clearly see that the recurrence is forward-stable (stable in increasing $p$) for positive $z$. For negative $z$, things are different. For $z<0$, $f_p$ as a function of the parameter $p$ for a fixed $z$ has a bump (it is not monotonous). See next graph.

Plotting the relative error gives:

One observes a loss of accuracy for increasing $p$. For $p$ above 27, you can see that the log plot of $f_p$ does not continue, simply because the recurrence gives negative values.
Step 2: can we do black-box numerical integration?
Again, for $z>0$, this would be meaningless since the recurrence works so well. We could try it for negative values of $z$. So let's have a look at the integrand:

Positive point: the integrand decreases very rapidly (after a peak). So one could approximate the infinite range integration to a finite range based on the value of the integrand.
Problematic point: the function is enormously peaked. Let's see if QUADPACK (used as back-end in scipy.integrate.quad
) can handle this. Note that I put the right point of integration arbitrarily at $r=200$. This could be optimized (if need be).
from scipy.integrate import quad
def F_neg(p,z):
Qres = np.zeros_like(z)
for i,zz in enumerate(z):
Q,err = quad(integrand,0,200,args=(float(p),float(zz)),epsrel=1e-12,epsabs=1e-12)
Qres[i] = Q
return Qres
Again, a plot of the relative error comparing this Python code to Maple results:

Perfect!
Preliminary conclusion: use the recurrence for $z>0$, use the numerical quadrature for $z<0$.
I don't think that using closed form solutions using hypergeometric functions (or LaguerreL functions) will help much since you will have linear combinations of these, possibly leading to cancellation problems.
Work to be done: check range of $z$ and $p$ for which the above holds. But that is for the poster of the question...