syms x
syms y
fx = x*exp(-0.5*(x+y)^2);
fy = y*exp(-0.5*(x+y)^2);
df = diff(fy,x)-diff(fx,y);
[x,y] = meshgrid(-3:0.1:3);
Df = eval(df);
surf(x,y,Df);