function Z=join(X,Y)
% join - Cartesian product of two arrays
%
% join :: [[N,M]], [[L,P]] -> [[N*L,M+P]].
%
% The result contains every pairing of rows from the two arguments

%Z=[repmat(X,size(Y,1),1) kron(Y,ones(size(X,1),1))];

% this is probably faster (no multiplications!)
m=size(X,1); n=size(Y,1);
Z=[repmat(X,n,1) Y(repmat(1:n,m,1),:)];
