Numerical approximation for derivative of matrix

I am trying to evaluate the shear stresses within a rectangular section, where the shear stress is given by the derivative of a variable phi. I divided the section by introducing 50 gridlines in the x and y axis, and got the following phi values:

In order to now get the shear stress, I need to apply the following formula:

where

However, I am unsure on what the best way is to get the two derivatives Txz and Tyz and combine them to get the shear force. Could someone lend some advice?