clear all;
m=7;n=12; % for example
% direct method giving U (with memory issues for big m and n):
B=round(9*rand(m,n))
D=kron(B,B);
U=D(:,n*(0:(n-1))+(1:n)); % straightforward "good column indices"
% indirect method giving V (without memory issues):
V=zeros(m^2,n);
for q=1:n
    V(:,q)=kron(B(1:m,q),B(:,q));% "1D" Kronecker products
end;
max(max(abs(U-V))),% checking U and V are identical