I want to differentiate a function `f(x, y)` with respect to `x` and `y` numerically using the difference quotient

```
def dfdx_v2(x, y, h):
    return 0.5 * (f(x+h, y) - f(x-h, y)) / h

def dfdy_v2(x, y, h):
    return 0.5 * (f(x, y+h) - f(x, y-h)) / h

```

where `h` is the length of the interval. In my assumption, the implementation above should be more and more accurate the smaller `h` gets. However, I noticed that this is not the case. When I plot the error for `h = [1e-12, 1e-4]` I get the following graph:

[![enter image description here](https://i.stack.imgur.com/smCRv.png)](https://i.stack.imgur.com/smCRv.png)

I would have expected, that the error gets smaller and smaller if I let `h` go to zero. However, for very small `h` the error increases substantially. Another interesting observation is that area where the error oscillates very strongly. Why does that happen? I would have expected, that the error is a monotonic increasing function without any fluctuations. Are there any ideas why the error behaves this way?

Here is the code that I used to produce the graph above:

```
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return x * x + y * x + x * y**2 + np.exp(x) + np.sin(x) * np.cos(x) + np.sin(y)

def dfdx_v1(x, y):
    return 2 * x + y + y**2 + np.exp(x) + np.cos(x)**2 - np.sin(x)**2

def dfdy_v1(x, y):
    return x + 2 * x * y + np.cos(y)

def dfdx_v2(x, y, h):
    return 0.5 * (f(x+h, y) - f(x-h, y)) / h

def dfdy_v2(x, y, h):
    return 0.5 * (f(x, y+h) - f(x, y-h)) / h

def main():
    x = 2.1415
    y = -1.1415

    h_min = 1e-12
    h_max = 1e-4
    n_steps = 10000

    h = np.linspace(h_min, h_max, n_steps)

    error = list()

    for i in range(n_steps):
        error.append(np.sqrt((dfdx_v1(x, y) - dfdx_v2(x, y, h[i]))**2 + (dfdx_v1(x, y) - dfdx_v2(x, y, h[i]))**2))

    plt.plot(h, error, linewidth=0.2, label="Error")
    plt.yscale("log")
    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

```

---