This is a cross-post from Math.SE since the question has received very little attention. If it is too trivial, or in other ways not suited for this site, please let me know and I will remove it.
Suppose $f : \mathbb{R}^n \rightarrow \mathbb{R}$ is a differentiable (for simplicity, let's assume ${C}^\infty$) function, and we would like to find the following partial derivatives numerically (at an arbitrary point):
$$ \frac{\partial^2 f}{\partial x_i \partial x_j}, \frac{\partial^3 f}{\partial^2 x_i \partial x_k}, \frac{\partial^3 f}{\partial x_i \partial x_j \partial x_k} $$
up to some order in the accuracy (say, $\mathcal{O}(\Delta x_i^{q_i}, \Delta x_j^{q_j})$ for the former, $\mathcal{O}(\Delta x_i^{q_i}, \Delta x_j^{q_j}, \Delta x_k^{q_k})$ for the latter, where $q_i, q_j, q_k \in \mathbb{N}$).
For $i = j = k$, I can just use the finite difference coefficients (see for instance the Wiki article), and then for $i \neq j \neq k$ I suppose I could iteratively apply the above provedure, but this may be numerically unstable.
I could instead simply write down the coefficients by hand (see for example here, here, or here for related questions), but this quickly becomes unfeasible and rather error-prone for higher orders in accuracy, especially for the third mixed derivative, so I'm wondering if there is a multi-parameter generalization (i.e. for $i \neq j \neq k$) to the above so I can algorithmically generate the coefficients to some given order in accuracy of each parameter.
For instance, if $n=2$, and I'm looking for $\partial_x \partial_y f$ up to $\mathcal{O}(\Delta x^2, \Delta y^2)$, I should get as output the 2-tuple of points $[(1, 1), (1, -1), (-1, 1), (-1, -1)]$ and their corresponding coefficients $[1/4, -1/4, -1/4, 1/4]$ (as per this question).
Any advice or references would be appreciated!
