subroutine leinit(n,xjac,wt,a,b)
*******************
*** Initialization for Legendre-Galerkin ******
***
***  Output:    a(k,j)=L_k(x_j)
***             B(k,j)=L_k(x_j)*w_j/(\gamma_k)
***             xjac:  Gauss-Lobatto points
***             wt:     weights for the discrete inner product
******************* 
      implicit double precision (a-h, o-z)
      dimension xjac(0:n),wt(0:n),a(0:n,0:n),b(0:n,0:n)
      alpha=0.d0
      beta =0.d0
      call jacobl(n,alpha,beta,xjac)

***   compute L_0(x_j) and L_1(x_j) ***         

      do i=0,n
         a(0,i)=1.d0
         a(1,i)=xjac(i)
      enddo

**** compute the rest of L_k(x_j) *****
      
      do i=0,n
         do j=1,n-1
            a(j+1,i)=((2*j+1.d0)*xjac(i)*a(j,i)-j*a(j-1,i))/(j+1.d0)
         enddo
      enddo
**** compute the weights ****

      a1=2.d0/(n*(n+1.d0))
      do i=0,n
         wt(i)=a1/(a(n,i)**2)
      enddo

**** compute B(k,j)=L_k(x_j)*w_j/(\gamma_k) ****

      do i=0,n-1
         do j=0,n
            b(i,j)=(i+.5d0)*a(i,j)*wt(j)
         enddo
      enddo
      do j=0,n
         b(n,j)=.5d0*n*a(n,j)*wt(j)
      enddo

      return
      end