Approximations of higher-order derivatives

Is your feature request related to a problem? Please describe.

I'm working on a project where it's important to estimate higher-order derivatives (e.g. 2nd, 3rd, 4th, and potentially mixed) of quantities in xarray datasets. xarray only has a helper for first derivatives (from https://github.com/pydata/xarray/pull/2398).

Describe the solution you'd like

I'd like to be able to call differentiate with a list of variables (e.g., data_array.differentiate(['x','x','y','y'])) and get the most accurate finite difference approximation of the derivative of the data array with respect to all of those variables.

Describe alternatives you've considered

It's tempting to chain together differentiate to approximate higher-order derivatives (e.g. data_array.differentiate('x').differentiate('x').differentiate('y').differentiate('y')). However, although I'm not an expert in finite difference approximations, I'm not sure if this is equivalent to computing the most accurate nth-order central (or average of forward and backward) differences, especially wrt. edges or with irregular grids. I could be wrong, though.

Additional context

The extent of my knowledge of how to estimate higher-order differences.

Start with the existing DataArray.differentiate API and the linked higher-order finite-difference reference. Determine how repeated and mixed variables should be represented, including irregular grids and edge behavior. Done means a defined approach can estimate the requested higher-order derivatives without relying on unresolved assumptions about chained calls.