clear; clf
syms x y;
Z=x.^2+sin(x).*y.^2
S=diff(Z,x)
[x,y]=meshgrid(-8:0.5:8)
mesh(x,y,eval(S)) 
%% 4-th: fmesh
clear; clf, syms x y; Z=x^2*y^3, 
Z1=diff(Z,x)
fmesh(Z1,[-1 2 -1 1]) %[ax bx ay by]
%% 5-th: mixed approach
clear; clf
syms x y;
Z=x^2+x*y^3
Z1=diff(Z,x)
[x,y]=meshgrid(-1:0.2:2,-1:0.2:2)
surf(x,y,eval(vectorize(char(Z1))))