Complex step differentiation is by now a well-known technique for numerically approximating derivatives of functions. Crudely put, one numerically approximates the derivative of a function evaluated at a real argument by taking a very tiny step in the imaginary direction. For instance, consider the following function (first studied by Lyness and Moler):
In[]:=
f[z_]:=Exp[z]/(Cos[z]^3+Sin[z]^3)
To approximate the derivative of this function at a real argument like π/5, one does the following:
With[{z=N[π/5],h=$MachineEpsilon},Im[f[z+h]]/h]
Out[]=
3.66101
The result compares favorably to the derivative computed symbolically:
N[f'[π/5]]
Out[]=
3.66101
Of course, to employ this method, one must be able to evaluate the function to be differentiated at complex arguments. Nevertheless, when one can do so, complex-step differentiation can give better results than the usual methods based on finite differences. More details on this method can be found in the paper of Squire and Trapp, and the paper of Martins, Sturdza, and Alonso.
However, if one wants to evaluate derivatives for complex arguments, one must resort to a different method. Very recently, a paper by Roelfs, Dudal, and Huybrechs shows that one can use a quaternionic step for evaluating complex derivatives. Although there are Mathematica implementations of quaternions in the standard packages and the Function Repository, it is more convenient to use the matrix representation of the quaternions to implement quaternionic step differentiation. In particular, we can use either the 2×2 or 4×4 matrix representations, along with MatrixFunction to evaluate the function at matrix arguments:
With[{z=N[π/4+π/3,20],h=10^-13},​​Tr[MatrixFunction[f,zIdentityMatrix[2]+h{{0,-1},{1,0}}].Inverse[{{0,-1},{1,0}}]]/(2h)]
Out[]=
3.1425957492811367264-2.8691752721699257125
With[{z=N[π/4+π/3,20],h=10^-13},​​Tr[MatrixFunction[f,zIdentityMatrix[4]+h{{0,0,-1,0},{0,0,0,1},{1,0,0,0},{0,-1,0,0}}].Inverse[{{0,0,-1,0},{0,0,0,1},{1,0,0,0},{0,-1,0,0}}]]/(4h)]
Out[]=
3.1425957492811367264-2.8691752721699257125
Compare with the symbolically-evaluated derivative:
N[f'[π/4+π/3],20]
Out[]=
3.1425957492811367264-2.8691752721699257125
The examples above used a particular unit quaternion (more specifically, the  unit quaternion) to approximate the derivative, but as the paper indicates, any unit quaternion of the form c1  + c2  + c3  (i.e., unit quaternions with zero real part) ought to do; that said, more investigation is needed on how to pick the best unit quaternion that gives an accurate derivative approximation, as well as whether a 4×4 might be needed when the 2×2 implementation does not suffice.

CITE THIS NOTEBOOK

Quaternionic step differentiation​
by J. M.​
Wolfram Community, STAFF PICKS, March 14, 2025
​https://community.wolfram.com/groups/-/m/t/3416627