For 1D this seems trivial. Something like x.^[0:N] should work.
Here's one way for 2D. It can be modified for 3D points or as high power as you want. This produces an inline function that takes an array of (x,y) points and outputs values of 2d monomials (up to power 4 in this particular code) at given points.
r = @(n) reshape(ones(n-1, n) + triu(ones(n-1, n), 1), [1, n-1, n]);
vander = @(x) [ones(size(x,1), 1) x ...
prod(reshape(x(:, r(3)), [size(x,1), 3, 2]), 3) ...
prod(reshape(x(:, r(4)), [size(x,1), 4, 3]), 3) ...
prod(reshape(x(:, r(5)), [size(x,1), 5, 4]), 3) ];
The lines above should be equivalent to
vander = @(x, y) [ones(size(x,1), 1) x y x.^2 x.*y y.^2 ...
x.^3 x.^2.*y x.*y.^2 y.^3 ...
x.^4 x.^3.*y (x.*y).^2 x.*y.^3 y.^4 ];
If you don't want it vectorized, or you use a for loop, then it can be simplified:
r = @(n) ones(n-1, n) + triu(ones(n-1, n), 1);
vander = @(x) [ones(size(x,1), 1) x ...
prod(x(:, r(3))) ...
prod(x(:, r(4))) ...
prod(x(:, r(5))) ];
I suppose you could also form 1D monomials for each dimension, tensor multiply them all together to form a N-D table and then take the triangular/tetrahedral shaped upper-left part of that table
