On the first part of your question, I think you should not call gradient, because convolution with the operators is what calculates the gradients. This makes more sense:
opx = [1,0,-1;0,0,0;1,0,-1]/4;
opy = [-1,0,-1;0,0,0;1,0,1]/4;
Gx = conv2(im, opx);
Gy = conv2(im, opy);
However, the choice of operators puzzles me. [1 0 -1; 2 0 -1; 1 0 -1] would be a more conventional choice for an x-gradient operator.
To get the elements of the Hessian, you can convolve Gx and Gy with gradient operators, something like
Gxx = conv2(Gx, opx);
and so on. Again, the motivation for the choice of operator isn't clear - the operators [-1 1] and [-1; 1] introduce a half-pixel offset into the results, so can be an awkward choice.
An alternative approach is just to use the gradient function instead of convolution. However, convolution gives you more control, and permits the introduction of some smoothing and the use of symmetrical operators.
Depending on your data, you may well need some smoothing as well as differencing to form estimates of gradients or Hessian elements.
