If you want a symbolic-like gradient you'll have to do it with symbolic variables:
syms x y
F = x^2 + 2*x*y − x*y^2
dF = gradient(F)
From there you might generate m-functions, see matlabFunction (If you don't have access to the symbolic toolbox look at the file exchange for a submission by John d'Errico that does symbolic operations (on polynomials only?))
If you want to work with numerical calculations in the end you can use the gradient function. But then I'd guess you'd need to vectorize your function:
f = inline(x.^2+3*x.*y,'x','y')
or preferably start using anonymous functions instead of inline:
f = @(x,y) x.^3 - 3*x.*y.^3
If you want to create a function for the gradient operator you might get something like this to do what you want:
gF = @(f,x,y) gradient(f(x,y),x(1,:),y(:,1))
But then you're venturing into a region where you should give the chebfun (also found on the file exchange) project some consideration.
HTH
