parameter (MM=2, NN=4)
       dimension al(0:NN), c(0:MM,0:NN,0:NN), Dm(0:NN,0:NN)

       pi=4.0*atan(1.0)
       do 5 j=0, NN
       al(j)=cos(pi*j/NN)
 5     continue
   
       do 10  j=0, NN
          zeta=al(j)
          c(0,0,0)=1
          c1=1

          do 20 n=1, NN
             c2=1
             do 30 nu=0, n-1
               c3=al(n)-al(nu)
               c2=c2*c3
               do 35 m=0, MM
                 if (m .gt. 0) then
                    ss=m*c(m-1,n-1,nu)
                 else
                    ss=0
                 endif
                 c(m,n,nu)=((al(n)-zeta)*c(m,n-1,nu)-ss)/c3
 35            continue
 30          continue          
             do 40 m=0, MM
                if (m .gt. 0) then
                   ss=m*c(m-1,n-1,n-1)
                else
                   ss=0
                endif
                c(m,n,n)=c1*(ss-(al(n-1)-zeta)*c(m,n-1,n-1))/c2
  40         continue
          c1=c2
  20      continue
       do 50 nu=0, NN
          Dm(j, nu)=c(MM,NN,nu)
  50   continue
  10   continue

       do 60 i=0,NN
        write (11, 99) Dm(i, 0), Dm(i,1), Dm(i,2), Dm(i,3), Dm(i,4)
 60    continue

 99    format (5(1x, f9.4))
       stop
       end