4th order accurate difference formula less accurate than 2nd order formula? I understand that there are different difference expressions for calculating numerical derivatives.
For example, this one is a central difference formula supposed to be 2nd order accurate, i.e. error decreases as O(h^2) 
$\frac {dy}{dx} = \frac {y(x+h)-y(x-h)}{2} $
This formula is supposed to be 4th order accurate, i.e. error decreases as O(h^4)
$\frac {dy}{dx} = \frac {y(x-2h)-8y(x-h)+8y(x+h)-y(x+2h)}{12} $
Is the second formula always better than the first formula or not ? I have heard someone say that if increment h is not small enough 2nd formula will be more inaccurate than first one. If this is the case, is there a way to figure out the limit of h below which 2nd formula will give more accurate values of derivatives compared to the 1st formula ? Thanks.
 A: Your formula seems of. Set
$$
Df(x,h)=\frac{f(x+h)-f(x-h)}{2h}
$$
then
$$
Df(x,h)=f'(x)+\frac{f'''(x)}{6}h^2+O(h^4)
$$
To eliminate the quadratic term, you have to compute
$$
4Df(x,h)-f(x,2h)=3f'(x)+O(h^4)
$$
which indeed results in
$$
\iff f'(x)=\frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h}+O(h^4)
$$
This is also called Richardson extrapolation (more commonly known in the context of Romberg integration). Naming the expression $D_2f(x,h)$, the next interpolation step would be
$$
16D_2f(x,h)-D_2f(x,2h)=15f'(x)+O(h^6),
$$
but one can also obtain that error order by combining the central differences for $h,2h,3h$.

Now to the errors: The function evaluation will always introduce some noise of the relative magnitude $\mu$, where $\mu$ is the machine constant of the floating point type. Thus your error is a combination of the approximation error and that noise, or about 
$$\frac{M_0\mu}h+M_3h^2$$ for the simple central difference which is smallest when both terms are about equal, which in most cases gives $h\sim\sqrt[3]\mu$ if the magnitudes $M_0$ of the function values and $M_3$ of the third derivative (and below $M_5$ of the fifth derivative) are about equal. For double this gives $h\sim 10^{-5}$ reducing the (relative) error to about $10^{-10}$ as a rather hard lower bound.
For the Richardson extrapolation formula, one gets in the same way the error about $$\frac{M_0\mu}h+M_5h^4,$$ which is minimal at about $h\sim\sqrt[5]\mu$, for double this is $h\sim 10^{-3}$, giving about $10^{-12}$ as the best relative error.
The following graph demonstrates the errors of the one-sided and central difference quotients and the first Richardson extrapolation.


To answer the question, at $h=10^{-5}$ where the simple central difference gives best results with error $10^{-10}$, the extrapolation formula will be dominated by the noise term $\mu/h$, which is also of size $10^{-10}$. The actual constants involved in the error terms are different, so that the second error might be slightly larger than the first.
