subroutine derivmatrix(nn,mm,x,dm,n2)
**********************************************************
*  Give a set of collocation points x(0:nn), we have the corresponding
*  Lagrange polynomials:  h_j(x), this routine compute
* 
*    dm(k,j,l)=d^l_x h_j(x)|_{x=x(k)} for l=1,..,mm
*
*   Input: nn, x(0:nn):  an arbitrary set of distinct collocation points 
*          mm: the highest derivative needed
*          n2: the first two dimensions of dm (n2 >=nn)
*
*  Output: dm: the derivative matrices asdescribed above
*
*    local   Working array: c
*********************************************************
      implicit double precision (a-h,o-z)
      dimension x(0:nn), Dm(0:n2,0:n2,mm), c(0:mm,0:nn,0:nn)
   
       do j=0, nn
          zeta=x(j)
          do m=1,mm
             c(m,0,0)=0.d0
          enddo
          c(0,0,0)=1.d0
          c1=1.d0

          do n=1, nn
             c2=1.d0
             do nu=0, n-1
               c3=x(n)-x(nu)
               c2=c2*c3
               do m=0, mm
                 if (m .gt. 0) then
                    ss=m*c(m-1,n-1,nu)
                 else
                    ss=0
                 endif
                 c(m,n,nu)=((x(n)-zeta)*c(m,n-1,nu)-ss)/c3
               enddo
            enddo
              do 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-(x(n-1)-zeta)*c(m,n-1,n-1))/c2
              enddo
          c1=c2
       enddo
       do m=1,mm
          do nu=0, nn
             Dm(j, nu,m)=c(m,nn,nu)
          enddo
       enddo

      enddo
      return
      end