# Backwards finite difference for mixed partials at higher order

I'm trying to understand what is the general method for calculating a backwards difference for a mixed partial of $n$ variables. Let's start with one variable:

The forward and backward finite differences are

$$f'(x) =\frac{f(x+h)-f(x)}{h}$$

and

$$f'(x) = \frac{f(x) - f(x-h)}{h}$$

However, when we go to mixed partials for $n=2$ variables discretized on a 2D grid with $h = x_{i+1}-x_{i}$, we can use (a central difference).

$$\left(\frac{\partial^2 u}{\partial x\partial y} \right)_{i,j} = \frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4 \Delta x \Delta y} + \mathscr{O}\left((\Delta x)^2, (\Delta y)^2 \right) \qquad (1)$$

How would one go about finding $\mathscr{O}\left((\Delta x)^2, (\Delta y)^2 \right)$? Furthermore, what about backwards difference? Notice the above central difference needs information of the $i+1,j+1$ point which is problematic near or approaching a singularity.

EDIT: Is this approached with a matrix method of the Taylor expansion of the mixed partials?

$$f(x+\Delta x, y+ \Delta y) = \sum_{i=0}^{\infty} \dfrac{(\Delta x)^i}{i!} \dfrac{\partial^{i}}{\partial x^i} f(x, y+ \Delta y) = \sum_{i=0}^{\infty} \dfrac{(\Delta x)^i}{i!} \sum_{j=0}^{\infty} \dfrac{(\Delta y)^j}{j!} \dfrac{\partial^{i+j}}{\partial y^j \partial x^i} f(x, y)$$