function ndf=cder(f,N,nd)
% multidomain partial derivative
% N order of Chebyshev polynomials (n of points inside subdomain is N+1)
% number of domains
% f input vector, function to be derived
% ndf output vector first derivative of function

D=get_cheby_matrix(N);  % get Chebyshev matrix

% Derivative Domain 1 
sdf=-D*f(1:N+1)';        % derivative in first domain
ndf(1:N+1)=sdf;

for id=1:nd-1,
    sf=f((id*(N-1)+1):(id*(N-1)+1+N));
    sdf=-D*sf(1:N+1)';        % derivative in subdomain
    ndf((id*(N-1)+3):(id*(N-1)+1+N))=sdf(3:N+1);
end