[X,Y] = meshgrid(0:.2:5);
Z1 = (Y.^2 + X.^2)/2;
Z2 = X.*Y;
[DX,DY] = gradient(Z2,.2,.2);
%contour(X,Y,Z1)
hold on
quiver(X,Y,DX,DY)
%colormap hsv

for i=-5:2:5
    if i<0
        x0=sqrt(-i)+.001;
    else
        x0=0;
    end
    x=linspace(x0,5,50);
    y=sqrt(x.^2+i);
    plot(x,y,'k')
end
axis([0,5,0,5])
axes square
hold off