Exa = 2;

switch Exa
case 1
  function y = f(x)
    y = exp(-x).*sin(x);
  endfunction
  a = 0; b = 1*pi;
  Aexact = 0.5 - 0.5*exp(-b)*(cos(b)+sin(b))
  x = linspace(a,b);
  A1 = trapz(x,f(x))
  [A1, A1-Aexact]
  A2 = cumtrapz(x,f(x));
  figure(1)
  plot(x,f(x),';function;',x,A2,';integral;')
  xlabel('x'); ylabel('Integral'); grid on
case 2
  A3 = quad(@f,a,b);
  [A3, A3-Aexact]
case 3
  function res = dHz(al,R,x,z)
    res = R*(R-x.*cos(al))/sqrt(R^2-2*x.*R.*cos(al)+x.^2 +z.^2).^3;
  endfunction
  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
  plot(x,Fz)
  xlabel('x'); ylabel('Fz'); grid on;
end