Your example actually only gave 1 component of the gradient, at an unspecified fixed value of x. The gradient of your scalar function f(x,y) is a vectorial quantity with components in both x and y directions, so at any given set of coordinates (x1,y1), there will be 2 gradient components w1 and v1.
There is a function "gradient" in base matlab since R2006a, if you want to compute your objective function on a grid (Matlab's definition of a well-defined grid) and compute the gradient numerically from that.
x = xMin:d:xMax
y = yMin:d:yMax
[X,Y] = meshgrid(x,y)
Z = f(X,Y)
[W,V] = gradient(Z,d)
This method depends on the resolution of the grid that you pre-compute on, as I believe the gradient function simply implements a centered finite difference on the grid, with maybe some special consideration at the edge.
If you are ultimately interested in the optimization of f(x,y), a non-gradient method for optimizing would probably be pretty easy/fast in 2 dimensions? There is "fminsearch" in base matlab that could be useful.
