What is different between gradient and imgradient functions in MATLAB
Question 1:
gradient returns both estimated horizontal $\left(\frac{\partial I}{\partial x}\right)$ and vertical components $\left(\frac{\partial I}{\partial y}\right)$ of the gradient vector, which in your example code are Gx and Gy. To do so, gradient computes for each point $X(i)$ the half difference of $X(i-1)$ and $X(i+1)$, which is the most basic gradient estimation method (it's called the central difference method). 
Then, computing

g = sqrt(Gx.^2+Gy.^2);  


gives you the gradient magnitude at each point.

imgradient differs since it directly outputs the gradient's magnitude (g) and direction (which you suppressed using ~) at each point. In addition, it computes the gradient using a two-dimensional method applying a sobel kernel to each point and its neighbors. Finally, this function allows you to specify which kernel to use to compute the numerical gradient if you don't want to use the Sobel estimator. You can for example choose to estimate the gradient with the central difference method as in gradient.

Question 2 :
Sobel method is just the simplest contour detection method providing good results for a low computation cost. But there are plenty of others, that are used depending on the application I guess.
Compared to the central difference method, it achieves better noise rejection using a different kernel with coefficients designed to do so. See as an example the kernel used for horizontal gradient estimation :
$$
\begin{matrix}
-1 &amp; 0 &amp; 1\\
-2 &amp; 0 &amp; 2\\
-1 &amp; 0 &amp; 1
\end{matrix}
$$
The considered point's row is emphasized (factor 2) compared to it's neighbors (factor 1). While the central difference method does not make that difference and weights all points with the same coefficient :
$$
\begin{matrix}
-1 &amp; 0 &amp; 1\\
-1 &amp; 0 &amp; 1\\
-1 &amp; 0 &amp; 1
\end{matrix}
$$