I think I figured out a solution. If there is a more elegant way, I am all ears. I want to extend this to five dimensions, and keeping track of the permutations is going to be tricky.
x=1:5;
y=1:6;
z=1:3;
X=repmat(x',[1 length(y) length(z)]);
Y=repmat(y,[length(x) 1 length(z)]);
Z=repmat(z,[length(x) 1 length(y)]);Z=permute(Z,[1 3 2]);
U=(X.^2+Y.^3.*Z).*exp(X.*Z.^2).*besselj(0,y);