function q = integralgl (fun, n)
%
% function q = integralgl (fun, n)
%
% Formula di quadratura di Gauss-Legendre in [-1,1] di ordine n <= 4.
if (n > 4)
  error('Formula di quadratura Gauss-Legendre non definita per n > 4');
end
xtable = [0, 0, 0, 0, 0;...
          -1/sqrt(3), 1/sqrt(3), 0, 0, 0;...
          -sqrt(3/5), 0, sqrt(3/5), 0, 0;...
          -sqrt(3/7+2/7*sqrt(6/5)), -sqrt(3/7-2/7*sqrt(6/5)),...
          sqrt(3/7-2/7*sqrt(6/5)), sqrt(3/7+2/7*sqrt(6/5)), 0;...
          -sqrt(5+2*sqrt(10/7))/3, -sqrt(5-2*sqrt(10/7))/3, 0,...
          sqrt(5-2*sqrt(10/7))/3, sqrt(5+2*sqrt(10/7))/3];
wtable = [2, 0, 0, 0, 0;...
          1, 1, 0, 0, 0;...
          5/9, 8/9, 5/9, 0, 0;...
          (18-sqrt(30))/36, (18+sqrt(30))/36,...
          (18+sqrt(30))/36, (18-sqrt(30))/36, 0;...
          (322-13*sqrt(70))/900, (322+13*sqrt(70))/900, 128/225,...
          (322+13*sqrt(70))/900, (322-13*sqrt(70))/900];
x = xtable(n + 1,1:n + 1)'
w = wtable(n + 1,1:n + 1)
q = w * fun(x);
