clf
z = -1:0.2:2; x = -0.4:0.1:0.7;
[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));   % define the anonymous function
  Hx(k) = quad(fx,0,2*pi)/(4*pi);  % integrate
  fz = @(al)dHz(al,1,x(k),z(k));   % define the anonymous function
  Hz(k) = quad(fz,0,2*pi)/(4*pi);  % integrate
endfor

subplot(1,2,1)
quiver(x,z,Hx,Hz,0.2)
grid on
axis('equal')


scal = 1./sqrt(Hx(:).^2+Hz(:).^2);  % length of each vector
Hx = scal.*Hx; Hz = scal.*Hz;       % normalize length

subplot(1,2,2)
quiver(x,z,Hx,Hz,0.1)
grid on
axis('equal')
