dHz = @(al,R,x,z)R*(R-x.*cos(al))./sqrt(R^2-2*x.*R.*cos(al)+x.^2 +z.^2).^3;
x = -0.5:0.05:0.8;
Fz = zeros(size(x));
for k = 1:length(x)
  fz = @(al)dHz(al,1,x(k),0);
  Fz(k) = quad(fz,0,2*pi)/(4*pi);
%%% pack it all in one line
%  Fz(k) = quad(@(al)dHz(al,1,0,z(k)),0,2*pi)/(4*pi);
end%for

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
end%for

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

