function dm=derivmatrix(nn,mm,x)
% **********************************************************
% *  Give a set of collocation points x(1:nn+1), 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(1:nn+1):  an arbitrary set of distinct collocation points 
% *          mm: the highest derivative needed
% *
% *  Output: dm: the derivative matrices asdescribed above
% *********************************************************

c=zeros(mm+1,nn+1,nn+1);

for j=1:nn+1
    zeta=x(j);
    c(1,1,1)=1;
    c1=1;
   
    for n=2:nn+1
        c2=1;
        for nu=1:n-1
            c3=x(n)-x(nu);
            c2=c2*c3;
            for m=1:mm+1
                if(m>1)
                    ss=(m-1)*c(m-1,n-1,nu);
                else
                    ss=0;
                end
                c(m,n,nu)=((x(n)-zeta)*c(m,n-1,nu)-ss)/c3;
            end
        end
        for m=1:mm+1
            if (m>1)
                ss=(m-1)*c(m-1,n-1,n-1);
            else
                ss=0;
            end
            c(m,n,n)=c1*(ss-(x(n-1)-zeta)*c(m,n-1,n-1))/c2;
        end
        c1=c2;
    end
    for m=1:mm
        for nu=1:nn+1
            dm(j,nu,m)=c(m+1,nn+1,nu);
        end
    end
end