The only difference I can see between your gradient3 function and MATLAB's gradient is that the latter returns the horizontal derivative as the first output, and your code returns as "x" derivative the vertical derivative.
Note that MATLAB arrays are stored such that the first index is vertical, and the second index is horizontal. Therefore, the code under ...
From MATLAB documentation, subsection Algorithms, the gradient G along the second axis of A is calculated as the central difference:
For j in 2:N-1, where N = size(A,2):
G(:,j) = 0.5*(A(:,j+1) - A(:,j-1));
At the edges where central difference cannot be used, the gradient is calculated by:
G(:,1) = A(:,2) - A(:,1);
G(:,N) = A(:,N) - A(:,N-1);
The gradient by definition is a vector formed by three directional derivatives, ∂F/∂x, ∂F/∂y and ∂F/∂z, so if you only need the z component, at least we have two perspectives, which are, in order of recommendation:
1.- You can calculate it directly with diff(F,z). Check the docs https://www.mathworks.com/help/symbolic/differentiation.html
2.- Or you can ...