% D=hefgsdiff(n,x) returns the first-order differentiation matrix of size
% n by n, at the Hermite-Gauss points x, which can be computed by 
% x=hegs(n), associated with Hermite function approach  
% Use the function: hefun() 
% Last modified on December 22, 2011

function D=hefgsdiff(n,x)

if n==0, D=[]; return; end;
  xx=x;y=hefun(n-1,xx); nx=size(x); 
if nx(2)>nx(1), y=y'; xx=x'; end; % xx, y are column-n vectors
  D=(xx./y)*y'-(1./y)*(xx.*y)';   % see (7.93) 
  D=D+eye(n);                     % add the identity matrix so that 1./D can be operated                                     
  D=1./D; D=D-eye(n);             % zero entry on the diagonal
   
 return;