Implementing the (polynomial) "kernel trick" in matlab

Given p = 2 vectors x in R^n, y in R^n, how to compute the kernel of order m > 2, which is of the form

A = [ones(n,1), x, y, x.*y, x.^2, y.^2, x.^2 *y, ..., x.^m, y.^m]  
q = nchoosek(p+m, m)

A in R^(nxq)

For example, given

x=[1 2]'
y=[2 1]'
m = 3;

Output:

A = [1 1 2 1 4 1 8;1 2 1 4 1 8 1]