[FX,FY] = gradient(F) returns the numerical gradient of the
matrix F. FX corresponds to dF/dx, the differences in x (horizontal)
direction. FY corresponds to dF/dy, the differences in y (vertical)
direction. The spacing between points in each direction is assumed to
be one. When F is a vector, DF = gradient(F) is the 1-D gradient.
So, since we're only using the first return value, we're telling surf to use dF/dx for the color. If we look at the values:
We see that dF/dx ranges from -0.0853 to 0.1922. Now, the surface is going to tell the axes what the range of the color data is. We can see that happen using the caxis function.
If there were other objects in the axes which had color data, then these values would get stretched to include all of them.
Next, the surface will get a colormap from the axes. We can can see that colormap with the colormap function.
That returns a 64x3 array of values. The columns are the red, green, and blue values of 64 different colors. The surface is going to map the value -0.0853 to the first of those colors, and the value 0.1922 to the last of those colors. Values in between these two will choose colors from the 64 using linear interpolation.
The colorbar function shows a visual representation of how the axes is telling the surface to transform colors. You can see 64 different colors (actually kind of hard to see because the default colormap is quite smooth) and the values they map to.
On the other hand, we've told the surface object to use the array Z for the Z values. These range from -0.4218 to 0.4218. You can see this range on the Z ruler on the left side of the axes.
This shows that the function we're using for Z can be different from the function we're using for color. We could put other functions into the color data. For example, if I wanted to see the 1st derivative of Z with respect to y instead of x, I would do it like this:
[dfdx,dfdy] = gradient(Z);
Notice how the high values of df/dy (the yellow colors) are in the places were I would be going uphill if I was walking in the positive Y direction. And the low values of df/dy (the blue colors) are in places where I would be going downhill.
Or I could show the magnitude of the gradient like this:
surf(x,y,Z,sqrt(dfdx.^2 + dfdy.^2))
In this case, the large values are in areas where both df/dx and df/dy are large. Also notice in this case, that the bottom of the colormap is at 0 instead of being negative. That's because sqrt(dfdx.^2+dfdy.^2) will always be positive.
Another thing that's important to note here is that the surface is using a single color value for each of those little quadrilaterals. We can see that if we use a smaller, less smooth colormap.
That's because the default value of the surface's FaceColor property is 'flat'. We have another option, which is 'interp'. This will interpolate the color values across the interior of the quadrilaterals, as you can see here.
surf(x,y,Z,sqrt(dfdx.^2 + dfdy.^2),'FaceColor','interp')
Does that help?