This function does exactly what Matlab KRON does, but for large full matrices, the engine uses BSXFUN to accelerate the calculation.
Another advantage is no intermediate large matrices are generated (four temporary arrays in case of KRON).
Here is the benchmark code and result:
mem = memory;
maxn = (mem.MaxPossibleArrayBytes/32)^0.25;
n = 10:10:maxn;
tic; K = kron(A,B); t1=min(t1,toc);
tic; K = kronecker(A,B); t2=min(t2,toc);
gain(end+1) = t1/t2;
fprintf('Size A/B Speed gain\n');
fprintf(' %02d %1.2f \n', [n; gain]);
Size A/B Speed gain