interp2 does not return derivative information at the interpolation points.
If you want to generate that at some general interpolation point, it would seem simplest to just use a finite difference approximation for the partial derivatives at that point. I would strongly suggest a central finite difference for accuracy, as long as spline interpolation is involved. If you were using linear interpolation, then just use a forward difference.
That is, for reasonably small dx, just evaluate the interpolant at two points.
df/dx ~ (f(x + dx, y) - f(x - dx, y))/(2*dx)
You can do the same for the derivative in y, holding x fixed at the current value..
df/dy ~ (f(x,y + dy) - f(x,y - dy))/(2*dy)
So to compute the gradient vector [df/dx,df/dy] just requires 4 function evaluations.
This is NOT an exact computation of the derivative at any point. It may even be moderately far off, if the function is highly nonlinear, but then any such estimate can be wild on a massively nonlinear surface.
If instead, you wanted to infer the derivatives at every node of the lattice used by interp2, then gradient will NOT provide the exactly correct derivatives as implicitly used by interp2. At best, it would return only an approximation. You could use the same solution I suggest above at each node of the lattice, as long as you used forward or central differences along the boundary of the lattice.
If you truly wanted the EXACT derivatives, you are probably out of luck, unless you were willing to code the spline equations up yourself. I suppose you could throw my derivest code at the surface to estimate the partial derivatives as well as possible. At least it tries to provide an error estimate too, to tell you if it thinks there is a problem.
[X,Y] = meshgrid(-1:.1:1);
Z = exp(X+2*Y);
So a very simple surface. Now lets estimate the gradient of the surface at (x,y) = (1.23,-0.35), fist using a central difference estimate.
delta = 1.e-8;
xy = [1.23,-0.35];
diff(interp2(X,Y,Z,xy(1)+ delta*[1 -1],xy(2),'spline'))/(2*delta)
diff(interp2(X,Y,Z,xy(1),xy(2)+ delta*[1 -1],'spline'))/(2*delta)
How would I do so using my derivest tools?
fun = @(xy) interp2(X,Y,Z,xy(1),xy(2),'spline');
[grad,err] = gradest(fun,xy)
So my gradest utility was a little more accurate, at least in the computation of the derivative of the splined surface. You cannot know the true derivative of a general unknown function without the functional form.