numerically approximating the 1st derivative with N points It is my understanding that if I have 3 equally spaced samples
$$f(x_0-\Delta x), f(x_0), f(x_0 + \Delta x)$$ 
then the best approximation of $\frac{df}{dx}$ is 
$$f_{\Delta3}(x_0) = \frac{f(x_0+\Delta x) - f(x_0-\Delta x)}{2\Delta x} $$
because this essentially ignores the 2nd derivative at $x_0$.
Does this generalize to N equally-spaced points?
I am particularly interested in $N=5$ but the "obvious" solution
$$f_{\Delta5}(x_0) = \frac{f(x_0+2\Delta x) + \frac{1}{2}f(x_0+\Delta x) - \frac{1}{2}f(x_0-\Delta x) - f(x_0+2\Delta x)}{5\Delta x}$$
doesn't seem to reduce the error from low-order polynomials, so I guess it is probably wrong.
For example if $f(x) = ax^3 + bx^2 + cx$ and we are looking for the derivative near $x_0=0$, then 
$$\begin{align}
f_{\Delta 3}(0) &= \frac{(a\Delta x^3+b\Delta x^2+c\Delta x) - (-a\Delta x^3+b\Delta x^2-c\Delta x)}{2\Delta x} \\
&= 2a\Delta x^2 + c
\end{align}$$
(the squared term drops out), but then
$$\begin{align}
f_{\Delta 5}(0) &= \frac{(8a\Delta x^3+4b\Delta x^2+2c\Delta x) + \frac{1}{2}(a\Delta x^3+b\Delta x^2+c\Delta x) - \frac{1}{2}(-a\Delta x^3+b\Delta x^2-c\Delta x) - ((-8a\Delta x^3+4b\Delta x^2-2c\Delta x))}{5\Delta x} \\
&= \frac{17}{5}a\Delta x^2 + c
\end{align}$$
and I would have expected the cubic term to drop out also with a higher-order approximation, if I got the weights correct, so it looks like I made a mistake.
Is there a way to fix?
 A: We know that the error of the central difference is an even expression in $Δx$,
$$
d_2(Δx)=\frac{f(x+Δx)-f(x-Δx)}{2Δx}=f'(x)+c_2(Δx)^2+c_4(Δx)^4+…
$$
Richardson extrapolation now tells you how to eliminate the first error term by combining the expressions for $Δx$ and $2Δx$, 
\begin{align}
d_4(Δx)&=\frac{4d_2(Δx)-d_2(2Δx)}3
\\[1.5em]
&=\frac43\frac{f(x+Δx)-f(x-Δx)}{2Δx}-\frac13\frac{f(x+2Δx)-f(x-2Δx)}{4Δx}
\\[1.5em]
&=\frac{-f(x+2Δx)+8f(x+Δx)-8f(x-Δx)+f(x-2Δx)}{12Δx}
\\[1.5em]
&= f'(x)+\tilde c_4(Δx)^4+…
\end{align}
A: To get the terms drop out in the most effective way we can proceed as follows.  We can formally write down the Taylor expansion of a function as:
$$\exp(h D)f(x) = f(x+h)$$
where $D$ is the differential operator. Then we can find many possible ways to express the derivative in terms of finite differences by writing the r.h.s. in terms of finite difference operators. E.g. we can write:
$$f(x+h) = (1+\Delta)f(x)$$
where $\Delta$ acts on $f$ as:
$$\Delta f(x) = f(x+h) - f(x)$$
We can thus write:
$$\exp(h D)f(x) = (1+\Delta)f(x)$$
This allows you to formally express the differential operator in terms of the finite difference operator. 
While this will yield a formally correct expression, it will be in terms of forward differences and will thus yield an asymmetric expression. What we want here is a symmetric expression, but this obtained in just the same say, you just consider the symmetric finite difference expression:
$$\left[\exp(h D)-\exp(-h D)\right]f(x) = f(x+h) - f(x-h)$$
We can write this as:
$$\sinh(h D)f(x) = \Delta_s f(x)$$
where $\Delta_s$ is the average between the forward and backward finite difference. We can thus formally write:
$$D = \frac{1}{h}\operatorname{arcsinh}{\left(\Delta_s\right)} = \frac{1}{h}\left[\Delta_s - \frac{\Delta_s^3}{6 } + \frac{3 \Delta_s^5}{40 }-\frac{5 \Delta_s^7}{112 }+\cdots\right]\tag{1} $$
To compute the powers of $\Delta _s$, it is convenient to introduce the shift operator $E$ that acts like:
$$E f(x) = f(x+h) $$
We can then write:
$$\Delta_s = \frac{E - E^{-1}}{2}$$
We can then compute the powers of $\Delta_s$ using the binomial expansion, e.g.:
$$\Delta_s^3 =\frac{1}{8}\left[E^3 - 3 E + 3 E^{-1} - E^{-3}\right]$$
If we do this in Eq. (1) to order $\Delta_s^7$, we obtain the expression:
$$D \approx \frac{1}{h}\left[\frac{1225}{2048}\left(E - E^{-1}\right) -\frac{245}{6144}\left(E^3 - E^{-3}\right)+\frac{49}{10240}\left(E^5 - E^{-5}\right) -\frac{5}{14336}\left(E^7 - E^{-7}\right)\right] $$
So, we have:
$$
\begin{split}f'(x)\approx \frac{1}{h}&\left[\frac{1225}{2048}\left(f(x+h) -f(x-h)\right) -\frac{245}{6144}\left(f(x+3 h) -f(x-3 h)\right)+\right.\\
&\left. \frac{49}{10240}\left(f(x+5h) - f(x-5 h)\right) -\frac{5}{14336}\left(f(x+7 h) -f(x-7 h)\right)\right]
\end{split}
$$
