I needed to find an O(h2) method to find f'''(x). Using Taylor expansions, I found:
$$f'''(x)=\frac{f(x+2h)-2f(x+h)-2f(x-h)+f(x-2h))}{2h^3} + O(h^2)$$.
However, I have also found that:
$$f'''(x)=\frac{f(x-3h)-6f(x-2h)+12f(x-h)-10f(x)+3f(x+h)}{2h^3} + O(h^2)$$ works too.
Are the formulas for approximation unique? Or can both of these be a good $O(h^2)$ approximation? If no, where have I made a mistake?
Also, I tried evaluating both of these formulas for $f(x)=\exp(x)$ with $h=10^-1, 10^-2,...,10^-9$ with $x=0$, but it produced huge values in the magnitude of $10^{27}$ and when I plotted them it didn't show anything. Reasoning through, it seems likely that the values would indeed be large when h is small as the denominator is $h^3$, but what does this mean for picking values of h? And surely there would still be a visible trend in the graph?
I then tried plotting the errors instead, as $f'''(0)=\exp(0)=1$, but these also did not plot nicely. I wanted to check that the formulas were indeed $O(h^2)$ by evaluating them both at $h$ and $\frac{h}{2}$ and seeing the difference in error, but this doesn't seem possible. What am I doing wrong?
More generally, for any function whose $5$'th derivative is bounded,
$ \sum_{j=1}^n c_j f(x + d_j h) = f'''(x) h^3 + O(h^5)$ if $$ \eqalign{\sum_{j=1}^n c_j &= 0\cr \sum_{j=1}^n d_j c_j &= 0\cr \sum_{j=1}^n d_j^2 c_j &= 0\cr \sum_{j=1}^n d_j^3 c_j &= 6\cr \sum_{j=1}^n d_j^4 c_j &= 0\cr}$$ For any $5$ distinct $d_j$, the Vandermonde matrix with entries $d_j^i$, $i=0..4$, $j=1..5$ is nonsingular, so there is a unique $(c_1, \ldots, c_5)$ that makes this work.
Your first formula is wrong. The second is correct, and is the only correct one that uses $f(x-3h), f(x-2h), f(x-h), f(x)$ and $f(x+h)$.
But you shouldn't try this numerically with extremely small $h$'s, because roundoff error becomes very bad.
EDIT: Here, for example, is a plot (done using Maple) of $\text{error}/h^2$ in the case $f = \exp$, $x = 0$. Roundoff error messes things up below about $h=0.005$.