The basic operation is to take half the difference between the two values on either side of the point you are considering. For example,
dx1(2,2) = 0.5 * (a(2,3) - a(2,1)) (i.e. 0.5*(9-7))
and
dx2(2,2) = 0.5 * (a(3,2) - a(1,2)) (i.e. 0.5*(7-3))
So for your example, dx1(1,2) = 0.5 * (2 - 1).
This has to be a little different for points on the edge of the matrix when a neighbouring value is not available. In these cases the single sided difference is taken, so for example
dx1(3,1) = a(3,2) - a(3,1) (i.e. 7-15)
If you are familiar with convolution, dx1 is just
conv2(a, [0.5 0 -0.5])
except for the left and right columns, and dy1 is
conv2(a, [0.5 0 -0.5].')
except for the top and bottom rows.
