x = -0.5:0.05:0.8;
Fz = zeros(size(x));
for k = 1:length(x)
  fz = @(al)dHz(al,1,x(k),0);     % define the anonymous function
  Fz(k) = quad(fz,0,2*pi)/(4*pi); % integrate
%%% pack it all in one line, not as readable, but possible
%  Fz(k)=quad(@(al)dHz(al,1,0,z(k)),0,2*pi)/(4*pi);
endfor

figure(1)
plot(x,Fz)
grid on
axis([-0.5 0.8 0 1.2]);
xlabel('position x'); ylabel('Field Fz');

x2= 1.2:0.05:3;
Fz2=zeros(size(x2));
for k=1:length(x2)
  Fz2(k)=quad(@(al)dHz(al,1,x2(k),0),0,2*pi)/(4*pi); % packed in one line
endfor

figure(2)
plot(x,Fz,x2,Fz2)
grid on
axis([-0.5 3 -0.6 1.2]);
xlabel('position x'); ylabel('Field Fz');

