%=====================================================
% Gauss - Hermite quadrature
%=====================================================
% n : number of points for evaluating f
%=====================================================

function [x,w] = gauss_her(n);


% Build Hn (Hermite Poly of order n) by recurrence :
h(1,1)=1;
h(2,1:2)=[2 0];
for k=2:n
    h(k+1,1:k+1)=2*[h(k,1:k) 0]-2*(k-1)*[0 0 h(k-1,1:k-1)];
end


Hn       = h(n+1,:);
Hn_deriv = polyder(Hn);

x        = roots(Hn);
w        = 2^(n+1)*factorial(n)*sqrt(pi)  ./...
           (polyval(Hn_deriv,x)).^2;