function [x,w,q]= qgaush(n,f,par)
%QGAUSH  Numerical evaluation of integrals (quadrature).
%        [x,w]=qgaush(n) computes the abscissas and weights 
%        of the Gauss-Hermite n-point quadrature formula with
%        integration limits (-inf,+inf).  
%        [x,w,q]=qgaush(n,f,par) computes the integral
%                  
%          q =   \int f(x) dx  ~=  \sum w[i] f(x[i],par)/exp(-x[i]^2)
%                                    i   
%        in addition to the abscissas and weights.  'f' is a 
%        string containing the name of the function.
%

%        Ellen McGrattan, 12-5-91
%        Revised, 2-22-95, ERM


% recursively compute coefficients of Hermite polynomials
p = 1;
q = [2;0];
for i=1:n-1;
  r = 2*[q;0]-2*i*[0;0;p];
  p = q;
  q = r;
end;

% the roots of the polynomials are the abscissas
x = sort(roots(r));

% derive the weights
d = [n:-1:0]'.*r;
d = d(1:length(d)-1);
p = 2;
q = [4;0]; 
for i=2:n-1;
  r = 2*[q;0]-2*i*[0;0;p];
  p = q;
  q = r;
end;
w = sqrt(pi)*polyval(r,x)./polyval(d,x);

% compute the integral if all input and output arguments are given
if nargin>1 & nargout>2;    
  x      = x(:);
  q      = w'*(feval(f,x,par)./exp(-x.*x));
end;

return

% an alternative way to get the weights:
p = 1;
q = [2;0];
A = [ones(1,n);
     2*x'];
for i=1:n-2;
 r = 2*[q;0]-2*i*[0;0;p];
 p = q;
 q = r;
 A = [A;polyval(r,x)']; 
end;
w = A\[sqrt(pi);zeros(n-1,1)]; 



