Mathworks is unfortunately often reluctant to disclose the details of their algorithms, and I'm afraid this seems to be the case for the 'divergence' function. I see no mention of the method used in the documentation for this function.
However, with a little effort you can probably determine the mathematical nature of their algorithm. For the purposes of making the call
div = divergence(X,Y,Z,U,V,W)
define X, Y, and Z with 'ndgrid' using monotone, non-uniformly spaced values for some reasonably small n by m by p size. Then first define U, V, and W linearly
U = A*X + B*Y + C*Z + D
V = E*X + F*Y + G*Z + H
W = I*X + J*Y + K*Z + L
for some arbitrary coefficients A, B, C, ..., L. Then 'div' will surely be equal to the constant sum A+F+K throughout the entire field, including at the edges. (If not, you can request your money back.)
Next change U, V, W to that of quadratic functions:
U = A*X.^2+B*Y.^2+C*Z.^2+D*Y.*Z+E*Z.*X+F*X.*Y+G*X+H*Y+I*J+K
and similarly with 20 more arbitrary coefficients in V and W. The exact divergence for these is of course
2*A*X+E*Z+F*Y+G + corresponding terms from V and W
and can be compared with matlab's 'div'. If they are equal (to within a reasonable round-off error,) that means Mathworks has used at least a second order approximation for their partial derivatives. If not, they are content with a first order approximation. If the answer was yes and you are ambitious, you can go on to higher order functions of X, Y, and Z for U, V, W and eventually pin down exactly what order of approximations they are using. Be sure to check the values on the faces, edges, and corners.
Of course, if you are lucky perhaps one of Mathworks' representatives might eventually decide to reveal at least that much information and save you all that trouble.
