I'm trying to make an algorithm that finds the first 10 or so terms of a function's Taylor series, which requires finding the nth derivative of the function for the nth term. It's easy to implement derivatives by following the definition of the derivative: $$f'(x) = \lim_{h\to0}\dfrac{f(x+h)-f(x)}{h}$$ implemented here in Python:

dx = 0.001
def derivative(f, x):
    return (f(x + dx) - f(x)) / dx

The value seems to be even closer to the actual value of the derivative if we define it like this:

dx = 0.001
def derivative(f, x):
    return (f(x + dx) - f(x - dx)) / (2 * dx)

which just returns the average of (f(x + dx) - f(x)) / dx and (f(x) - f(x - dx)) / dx.

For higher order derivatives, I implemented a simple recursive function:

dx = 0.001
def nthDerivative(f, n, x):
    if n == 0:
        return f(x)
    return (derivative(f, n - 1, x + dx) - derivative(f, n - 1, x - dx)) / (2 * dx)

I tested the higher order derivatives of $f$ at $1$, where $f(x)=x^9$, and as can be proved by induction, $$\dfrac{d^n}{dx^n}(x^k)=\dfrac{k!}{(k - n)!}x^{k-n}$$

Therefore, the nth derivative of $f$ at $1$ is $\dfrac{9!}{(9 - n)!}$. Here are the values returned by the function for n ranging from 0 to 9:

n             Value  Intended value
-----------------------------------
0             1.000               1
1             9.000               9
2            72.001              72
3           504.008             504
4          3024.040            3024
5         15120.252           15120
6         60437.602           60480
7         82298.612          181440
8      32278187.177          362880
9   95496943657.736          362880

As you can see, the values are waaaay off for $n$ greater than $5$.

What can I do to get closer to the actual values? And is there an algorithm for this that doesn't have $O(2^n)$ performance like mine?