Do you refer to the the Jacobian matrix of a function f:R^n->R^m? In this case you can program this code for an arbitrary function:
J[fct_List, var_List] :=
Module[{m, n}, m = Length[fct]; n = Length[var];
Simplify[Table[D[fct[[i]], var[[j]]], {i, 1, m}, {j, 1, n}]]]
In the case of your function a(x,y,z) you define it as follows:
r:={x, y, z};
a:=-my r/(r.r)^(3/2);
An then you simply enter
J[a,r]
