Approximating a function's derivative in machine precision Let's say I have a very complex analytical function $f$, in which finding its derivative at point $Q$, i.e. $f'(Q)$ is not practical. So we resort to finding the derivative numerically without analytically differentiate the function at point $Q$ (note that $Q$ is only a point, not an array).
Method 1: I have tried the $d$ order of finite difference. This doesn't give me a machine precision result when the order is too high or singularity behaviour occurs when steps are too small.
Method 2: Interpolation method. This doesn't give machine precision accuracy neither when differentiate the respected interpolant at point $Q$. i.e. $f'(Q) ≅ \sum_i L_i'(Q)f(x_i)$ where $L_i(Q)$ is the interpolant with basis $i$ and $x_i$ is the nodal point.
Method 3: Spectral method such as Legendre expansion and the like. No machine precision.
Method 4: least square polynomial fit. Just no.
Method 5: minimax polynomial approximation. Seems alright with around 13 decimal places accuracy in matlab.
Any better method? Suggestion? Alternatives?
I will be much appreciated. Thank you all!

If the precision is the priority (and speed is not) then consider what follows.
To increase the precision of numerical differentiation do the following:
1) Chose your favorite high-precision "standard" method based on some step size H.
2) Compute the value of the derivative with the method chosen in 1) many times with different step sizes h. Each time you may pick h as a random number from the interval (0.5*H/10, 1.5*H/10) where H is an appropriate step size for the method you use (yes, you should use smaller step sizes then the "optimal" one).
3) Average the results.
Your result may gain 2-3 orders of magnitude in the absolute error wrt. the non-averaged result.
https://arxiv.org/abs/1706.10219
