clf
x = -0.3:0.05:0.3; z=x;
h = 0.5; % optimal distance for Helmholtz configuration

[xx,zz] = meshgrid(x,z);
x = xx(:); z = zz(:);
Hx = zeros(size(x)); Hz = Hx;

for k = 1:length(x)
  fx = @(al)(dHx(al,1,x(k),z(k)-h)+dHx(al,1,x(k),z(k)+h));
  Hx(k) = quad(fx,0,2*pi)/(4*pi);
  fz = @(al)(dHz(al,1,x(k),z(k)-h)+dHz(al,1,x(k),z(k)+h));
  Hz(k) = quad(fz,0,2*pi)/(4*pi);
endfor

quiver(x,z,Hx,Hz,0.001)
grid on
