% example illustrating the divergence operator
clear, clf;
alpha=3;
[x,y] =meshgrid( -2 : .5 : 2, -2 : .5 : 2 );
M=length(x); N=length(y);
r2=x.^2+y.^2;
for k=1:M
  for l=1:N
    Ax(k,l)=exp(-r2(k,l)/alpha.^2)*x(k,l);
    Ay(k,l)=exp(-r2(k,l)/alpha.^2)*y(k,l);
  end
end
quiver(x,y,Ax,Ay)
xlabel('$$x$$','Interpreter','Latex','Fontsize',14)
ylabel('$$y$$','Interpreter','Latex','Fontsize',14)
title('$${\bf{f}}(x,y) = xe^{-(x^2+y^2)/9}\hat{\bf{x}}+ye^{-(x^2+y^2)/9}\hat{\bf{y}}$$','Interpreter','Latex','Fontsize',14)
figure;
D=divergence(x,y,Ax,Ay);
surf(x,y,D)
xlabel('$$x$$','Interpreter','Latex','Fontsize',14)
ylabel('$$y$$','Interpreter','Latex','Fontsize',14)
zlabel('div$${\bf{f}}(x,y)$$','Interpreter','Latex','Fontsize',14)
title('$${\bf{f}}(x,y) = xe^{-(x^2+y^2)/9}\hat{\bf{x}}+ye^{-(x^2+y^2)/9}\hat{\bf{y}}$$','Interpreter','Latex','Fontsize',14)