Gypaets/trigradient​2

[ZX, ZY, ZXX, ZYY, ZXY] = trigradient2(X, Y, Z, T, M)
The derivatives of the function Z(X,Y) are calculated with a least squares linear regression. The system of equations is set up with Taylor series from each point to the adjacent vertices. If a vertex is connected to less than five vertices then vertices in two-edge distance are used too.
This derivation method provides much better results than first order approaches. Specially the error of the calculated second order field derivatives is significantly smaller than deriving two times the field with a first order function.
Input:
X= Vector with x-coordinates.
Y= Vector with y-coordinates.
Z= Matrix with function values on each point. If Z has multiple
columns the derivatives for each column are calculated.
Optional argument:
T= Triangulation (Nx3 matrix with polygon vertices). If not given the delaunay triangulation of X,Y is used.
M= Method used for calculation. Default value is 0.
0: One large equation system. Fast.
1: Multiple small equations systems. Slower but depending on the input values more accurate.

Output:
ZX=dz/dx
ZY=dz/dy
ZXX=d^2z/dx^2
ZYY=d^2z/dy^2
ZXY=d^2z/(dx dy)