subroutine orthzeroweights(n,win,a,b,c,z,w)
*******************************************
* Compute zeros of the orthogonal polynomial p_n defined by the
* three-term recurrence relation:
*      p_{k+1}(x)=(a(k) x -b(k)) p_k(x) -c(k) p_{k-1}(x), k>=0
*
*   Input: n: degree of polyomial
*          a,b,c: coefficients in the three term relations
*          win=\int w(x) dx (where w(x) is the weight functions)
*          
* Output: z, the n zeros of p_n.
*         w, the corresponding weights of the quadrature formula
********************************************
      implicit double precision (a-h,o-z)
      dimension a(0:n-1),b(0:n-1),c(0:n-1),z(n),w(n),eig(n,n),wk(2*n)

*** compute alp_j =b_j/a_j, j>=0;  beta_j=sqrt(c_j/(a_{j-1} a_j)), j>=1

      do i=0,n-1
         z(i+1)=b(i)/a(i)
      enddo

      do i=1,n-1
         b(i-1)=sqrt(c(i)/(a(i-1)*a(i)))
      enddo
      
      call dstev('v',n,z,b,eig,n,wk,info)

**** compute the weights
      do i=1,n
         w(i)=win*eig(1,i)**2
      enddo
      return
      end