function dp=mpolydif(p,i)
% differentiate multivariate polynomial wrt variable i
% see mpolyval for representation of polynomials
n=size(p);
m=length(n);
p3=reshape(p,prod(n(1:i-1)),n(i),prod(n(i+1:m)));
% now the 2nd dimension of p3 is index i