Numerical method to solve $\sum\limits_{k=0}^{n}\frac{1}{2k+1} \lim\limits_{\sigma \to 0} \frac{d^{2k}}{d\sigma^{2k}} \left(f(\sigma,\alpha)\right)$ As part of a larger physical model I am currently searching for a solution to the following expression, a numerical solution is fine as I am ultimately really after the numerical result. $\alpha$ is a physical value and I need solutions for upto $n = 200000$
$$ \sum\limits_{k=0}^{n}\frac{1}{2k+1} \lim_{\sigma \to 0} \frac{d^{2k}}{d\sigma^{2k}} \left(\exp(0.5[\alpha-\sigma]^2) . \mathrm{erfc}([\alpha-\sigma]/\sqrt2) \right), \alpha \in R $$
$\mathrm{erfc}$ is the complementary error function.
Obviously an iterative solution would be best as it allows to find solutions for all $n$ in one run.
 A: The function $f(\alpha-\sigma) = \exp\left(\frac{1}{2}\left(\alpha-\sigma\right)^2\right) \operatorname{erfc}\left(\frac{\alpha-\sigma}{\sqrt{2}}\right)$ is a particular case of Hermite function $H_\nu\left(x\right)$, namely $f(x) = H_{-1}\left(\frac{x}{\sqrt{2}}\right)$. It's symbolic order derivative is known in closed form:
$$
     \frac{d^n}{dx^n} H_\nu\left(x\right) = (-2)^n (-\nu)_n H_{\nu-n}\left(x\right) 
$$
where $(a)_n$ is the Pochhammer symbol. Thus
$$
    \lim_{\sigma\to 0}\frac{d^{2k}}{d \sigma^{2k}} f(\alpha-\sigma) = \frac{d^{2k}}{d \alpha^{2k}} f(\alpha) = 2^k (2k)! H_{-2k-1}\left(\frac{\alpha}{\sqrt{2}}\right)
$$
The expression you are after, i.e.
$$
   S_n = \sum_{k=0}^n \frac{2^k (2k)!}{2 k+1} H_{-2k-1}\left(\frac{\alpha}{\sqrt{2}}\right)
$$
satisfies a recurrence equation (found with Mathematica):
$$
    \left(5 \alpha ^2+2 n \left(\alpha ^2+4 n (n+6)+52\right)+81\right) S_{n+1}+\left(-5
   \alpha ^2-2 n \left(\alpha ^2+4 n+20\right)-52\right) S_{n+2}-2 (n+2) (2 n+3)^2
   S_n+(2 n+7) S_{n+3}=0
$$
with self-evident initial conditions. The recurrence equation can be a basis for an efficient implementation:

