subroutine heinit(n,x,a,grw)
******************************************
* Initialization for Scaled Hemite functions U_i(x)=alpha_i H_i e^(-x*x/2)
* such that \int U_i u_j dx=delta_ij
* 
* Input: x(j), j=1,n: Hermit-Gauss points
* Output: a(j,i)=U_i(x_j):=alpha_i H_i(x(j))*exp(-x(j)*x(j)/2)
*         grw(i)= scaled Hermite-Gauss weight
*
*  It works for n up to at least 1400 ! It will eventually overflow for n very
*  large...
*  To make it work for even larger n, one should do it as described in
*  my SIAM laguerre paper, see ~/laguerre/lib/lainit.f !
*****************************************
      implicit double precision (a-h,o-z)
      dimension x(n),a(n,0:n-1),grw(n)
      spi=sqrt(acos(-1.d0))

      do j=1,n
         t0=1.d0/sqrt(spi)
         a(j,0)=t0*exp(-.5d0*x(j)*x(j)/2)
         a(j,1)=2*x(j)*a(j,0)/sqrt(2.d0)
         
         do i=2,n-1
             a(j,i)=2*x(j)*a(j,i-1)/sqrt(2.d0*i)
     1             -2*(i-1)*a(j,i-2)/sqrt(4.d0*i*(i-1))
          enddo
          do i=0,n-1
             a(j,i)=a(j,i)*exp(-.5d0*x(j)*x(j)/2)
          enddo
      enddo

**** Compute the scaled Gauss-Radau weight

      do j=1,n
         grw(j)=1.d0/(n*a(j,n-1)**2)
      enddo
      return
      end