[X,Y] = meshgrid(-2.2:.01:2.2);
%Z = X.*exp(-X.^2 - Y.^2);
Z=1./sqrt(X.^2+Y.^2);
% Z=0.5./sqrt((X-1).^2+Y.^2)+0.5./sqrt((X+1).^2+Y.^2);
[DX,DY] = gradient(Z,.01,.01);
DDY=-DY(1:20:end,1:20:end);
DDX=-DX(1:20:end,1:20:end);
XX=X(1:20:end,1:20:end);
YY=Y(1:20:end,1:20:end);

DDX(11:13,12)=0;
DDX(12,11:13)=0;
DDY(11:13,12)=0;
DDY(12,11:13)=0;
% 
% DDY=DDY(1:2:end,1:2:end);
% DDX=DDX(1:2:end,1:2:end);
% XX=XX(1:2:end,1:2:end);
% YY=YY(1:2:end,1:2:end);

figure
% contour(X,Y,Z)
quiver(XX,YY,DDX,DDY)
axis([-2.2 2.2 -2.2 2.2])
xlabel('x')
ylabel('y')
pause
hold on
contour3(X,Y,Z,[0.5 0.6 0.8 1.2 2 3.5])
zlabel('“dˆÊ')
grid
axis([-2.2 2.2 -2.2 2.2 0 3.5])
%quiver(X,Y,DX,DY)
hold off