% M3eP381.m
[x,y] = meshgrid(-2:0.25:2);
f = x.*exp(-((x-y.^2).^2+y.^2));
dx = x(1,2) -x(1,1); dy = y(2,1) - y(1,1);
[df_dx, df_dy] = gradient(f, dx, dy);
subplot(2,1,1)
contour(x,y,f), xlabel('x'), ylabel('y'), ...
hold on, quiver(x,y,df_dx, df_dy), hold off
subplot(2,1,2)
mesh(x,y,f), xlabel('x'), ylabel('y'),zlabel('f')