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_{hto0}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?