Numerical partial derivatives without using gradient

I am trying to implement numerical partial differentiation in MATALB wihtout using the built-in functions. I have a function vel = x+exp(-((x-x(1)).^2+(y-y(1)).^2)) of two variables with x ranging from -1 to 1 and y range from -2 to 2, with increments of 0.1 on both.

This is my code so far, but I am not sure it is correct:

x = [-1:0.1:1];
y = [-2:0.1:2]'; % Is this a legal operation? (the transpose) 

[X,Y] = meshgrid(x,y);
vel = X+exp(-((X-x(1)).^2+(Y-y(1)).^2)); % My function

derivx = (vel(2:end,1:end-1) - vel(1:end-1,1:end-1))./(x(2:end)-x(1:end-1));
derivy = (vel(1:end-1,2:end) - vel(1:end-1,1:end-1))./(y(2:end)-y(1:end-1));

It seems that when I check it with WolframAlpha derivx looks like the partial with respect to y and derivy looks like the partial with respect to x? Where is the mistake?

How would I go about to do the second order partial derivatives?

Thanks in advance for your help!

[Bjorn Gustavsson]

Yes, you are correct - and welcome to matlab...

For your purposes matlab stores matrices as

M(y,x)

So you are correct regards to the X-Y derivatives, and the only thing you have to do is to adjust for the row-column ordering. It took me a good while to get used to this and sometimes I still slip up...

HTH