Question
Let $X$ be a random variable for which we only have the value of its Moment Generating Function $M_X$ on a discrete set of points, I am looking for a stable method to compute: $$\frac{M_X(s) - M_X(t)}{s-t}=\frac{\mathbb{E}[e^{sX}-e^{tX}]}{s-t}, \qquad s \approx t \approx 0.$$
Thoughts/ Attempts
The problem is, clearly, that both the numerator as the denumerator becomes small as $|s-t| \rightarrow 0$. We can however use the Taylor expansion of $e^{sX}$ to overcome this problem, indeed we have (with $ \gamma_n(x,y) = \sum_{m=0}^{n-1} x^m y^{n-1-m} $): \begin{align*} \frac{\mathbb{E}[e^{sX}-e^{tX}]}{s-t} &= - \sum_{n=1}^\infty \frac{\gamma_{n-1}(s, t)}{n!} \left( \frac{d^n M_X}{ds^n} \right)\bigg|_{s=0}. \end{align*} This resolves our original problem. However it introduces another problem, namely the computation of: $\left( \frac{d^n M_X}{ds^n} \right)\bigg|_{s=0}$, which again causes numerical instability. I have tried to stabilize this differentiation. I see no method to do this (except by applying Cauchy's integral formula, which can not be applied here as we only know the value of $M_X(s)$ for real values $s$). Maybe we can rewrite this formula again in function of $M_X(s)$ but I am not sure how.
I have also implemented the suggestions found here but high order derivatives unfortunately can not be computed this way.
Background/Reason For Question
I am using this quantity in a recursion and numerical errors introduced by it make the recurrence fail. This happens after $10-20$ steps, some precision is lost in each subsequent step and this inaccuracy explodes..