Multinomial Expansion

MULTINOMIAL_EXPAND determines the matrix of powers for a multinomial expansion
of the form (x_1 + x_2 + x_3 + ... + x_ndim)^pow
Nmatrix - matrix of powers, each row representing a single term in the expansion
for example, the row [0,1,0,2] would represent (x_2)*(x_4)^2
Note, this is equivalent to finding all multiindices
k = [k_1,k_2,...,k_ndim] with |k|=sum(k)=pow
This can be useful for sparse-grid integration methods that use multiindices.

Ncoef - vector of coefficients (multinomial coefficient)

Ex. Evaluating the multinomial at a point x = [x_1, x_2, x_3, ... ]
-> sum(Ncoef .* repmat(x,size(Nmatrix,1),1).^Nmatrix)

Ex. Compute all multiindices of length 6 and order 4:
-> Nmatrix = multinomial_expand(4,6);

The method is recursive, so it can be a bit slow for smaller problems, but
it can work with relatively large size inputs (e.g. ndim = 50,100,...)
compared to previous versions (e.g. another implementation could
not handle ndim>15 due to memory constraints using a huge kron).
The savings in memory is because we only compute the
needed number of terms and don't use kron or factorial.