Hi Kevin,
It turns out that nDfy only gives the correct result because of the symmetry.
xo = 0.8;
yo = 0.2;
h = 0.01; %grid size deltax = deltay = 0.01
[nx,ny] = meshgrid(0:h:2); %x and y range from 0 to 2, step size h
z = nx.*exp(-(nx.^2+ny.^2));
%locate index values for xo = 0.8 and yo = 0.2.
nxo = find(nx(1,:)==xo);
nyo = find(ny(:,1)==yo);
[nDfx, nDfy] = gradient(z,h); %calculate the partials w.r.t x and y
nDfx(nxo,nyo) %incorrect value
nDfx(nyo,nxo) %correct value, but why?
nDfy(nxo,nyo) % happens to be the correct value
nDfy(nyo,nxo) % corret value
When we use meshgrid, the x-values vary across the 2nd dimension (horizontal) and the y-values vary down the 1st dimension (vertical). Sounds counterintutive, but that's how it works. Similarly, gradient assumes the x-direction is horizontal across the array and the y-direction is vertical down the array, which is consistent with the meshgrid convention. But that means that indexing into the results also has to reverse as well, which is what we're seeing.
Consider a simpler problem
x = 1:3;
y = 101:103;
[X,Y] = meshgrid(x,y);
Z = 5*X + Y
If we want the value of Z that corresponds to the second element of x and the first element y (z = 5*2 + 101 = 111), we'd need Z(1,2)
