CJam, 13 bytes

{ffff*::.+:~} 

This is an unnamed block that expects two matrices on top of the stack and leaves their Kronecker product in their place.

Test suite.

Explanation

This is just the Kronecker product part from the previous answer, therefore I'm here just reproducing the relevant parts of the previous explanation:

Here is a quick overview of CJam's infix operators for list manipulation:

• f expects a list and something else on the stack and maps the following binary operator over the list, passing in the other element as the second argument. E.g. [1 2 3] 2 f* and 2 [1 2 3] f* both give [2 4 6]. If both elements are lists, the first one is mapped over and the second one is used to curry the binary operator.
• : has two uses: if the operator following it is unary, this is a simple map. E.g. [1 0 -1 4 -3] :z is [1 0 1 4 3], where z gets the modulus of a number. If the operator following it is binary, this will fold the operator instead. E.g. [1 2 3 4] :+ is 10.
• . vectorises a binary operator. It expects two lists as arguments and applies the operator to corresponding pairs. E.g. [1 2 3] [5 7 11] .* gives [5 14 33].

ffff* e# This is the important step for the Kronecker product (but e# not the whole story). It's an operator which takes two matrices e# and replaces each cell of the first matrix with the second matrix e# multiplied by that cell (so yeah, we'll end up with a 4D list of e# matrices nested inside a matrix). e# Now the ffff* is essentially a 4D version of the standard ff* idiom e# for outer products. For an explanation of ff*, see the answer to e# to the Kronecker sum challenge. e# The first ff maps over the cells of the first matrix, passing in the e# second matrix as an additional argument. The second ff then maps over e# the second matrix, passing in the cell from the outer map. We e# multiply them with *. e# Just to recap, we've essentially got the Kronecker product on the e# stack now, but it's still a 4D list not a 2D list. e# The four dimensions are: e# 1. Columns of the outer matrix. e# 2. Rows of the outer matrix. e# 3. Columns of the submatrices. e# 4. Rows of the submatrices. e# We need to unravel that into a plain 2D matrix. ::.+ e# This joins the rows of submatrices across columns of the outer matrix. e# It might be easiest to read this from the right: e# + Takes two rows and concatenates them. e# .+ Takes two matrices and concatenates corresponding rows. e# :.+ Takes a list of matrices and folds .+ over them, thereby e# concatenating the corresponding rows of all matrices. e# ::.+ Maps this fold operation over the rows of the outer matrix. e# We're almost done now, we just need to flatten the outer-most level e# in order to get rid of the distinction of rows of the outer matrix. :~ e# We do this by mapping ~ over those rows, which simply unwraps them. 

MATLAB / Octave, 83 42 Bytes

Saved 41 bytes, thanks to FryAmTheEggman!

@(A,B)cell2mat(arrayfun(@(n)n*B,A,'un',0)) 

Test it here!

Breakdown

arrayfun is a disguised for-loop that multiplies n*B, for a variable n defined by the second argument. This works because looping through a 2D matrix is the same as looping through a vector. I.e. for x = A is the same as for x = A(:).

'un',0 is equivalent to the more verbose 'UniformOutput', False, and specifies that the output contains cells instead of scalars.

cell2mat is used to convert the cells back to a numeric matrix, which is then outputted.

Pyth, 14 12 11 bytes

JEsMs*RRRRJ 

Translation of Jelly answer, which is based on Büttner's Algorithm (ü pronounced when trying to make an ee sound [as in meet] in the mouth-shape of an oo sound [as in boot]).

Try it online (test case 1)!

Bonus: calculate B⊗A in the same number of bytes

JEsMs*LRLRJ 

Try it online (test case 1)!

Julia, 40 39 37 bytes

A%B=hvcat(sum(A^0),map(a->a*B,A')...) 

Try it online!

How it works

• For matrices A and B, map(a->a*B,A') computes the Kronecker product A⊗B.

  The result is a vector of matrix blocks with the dimensions of B.

  We have to transpose A (with ') since matrices are stored in column-major order.

• sum(A^0) computes the sum of all entries of the identity matrix of A's dimensions. For an n×n matrix A, this yields n.

• With first argument n, hvcat concatenates n matrix blocks horizontally, and the resulting (larger) blocks vertically.

APL (Dyalog Unicode), 10 bytes

{⍪/,/⍺×⊂⍵} 

Try it online!

How it Works

{ ... } ⍝ dfn that takes ⍺=A and ⍵=B ⍺×⊂⍵ ⍝ Product of each element of A with all of matrix B ⍝ Gives a nested array: an matrix of matrices ./ ⍝ Join rows ⍪/ ⍝ Join columns 

R, 98 86 bytes

function(x,y,`+`=array)aperm(apply(x,1:2,`*`,y)+c(w<-dim(y),v<-dim(x)),c(1,3,2,4))+v*w 

Try it online!

Reimplementation of .kronecker and outer for matrices. I do think there's a golfier approach out there, maybe using apply? 6 bytes golfed using apply and array thanks to Dominic van Essen!

The builtins are %x% for kronecker(A,B,"*") and %o% for outer(A,B,"*").

R, 120 bytes

function(A,B){l=list() a=dim(A) for(i in 1:a[2]-1)l[[i+1]]<-do.call(rbind,lapply(A,"*",B)[1:a+a[2]*i]) do.call(cbind,l)} 

Try it online!

Naive approach: calculate the subarrays \$a_{ij}B\$ and bind them together in the appropriate order. There are probably golfs here too, but I don't think it'll be shorter than the one above.

J, 10 bytes

This is one possible implementation.

[:,./^:2*/ 

J, 13 bytes

This is a similar implementation, but instead uses J's ability to define ranks. It applies * between each element on the LHS with the entire RHS.

[:,./^:2*"0 _ 

Usage

 f =: <either definition> (2 2 $ 1 2 3 4) f (2 2 $ 5 6 7 8) 5 6 10 12 7 8 14 16 15 18 20 24 21 24 28 32 (2 1 $ 1 2) f (1 2 $ 1 2) 1 2 2 4 2 f 5 10 

JavaScript (ES6), 79

Straightforward implementation with nested looping

(a,b)=>a.map(a=>b.map(b=>a.map(y=>b.map(x=>r.push(y*x)),t.push(r=[]))),t=[])&&t 

Test

f=(a,b)=>a.map(a=>b.map(b=>a.map(y=>b.map(x=>r.push(y*x)),t.push(r=[]))),t=[])&&t console.log=x=>O.textContent+=x+'\n' function show(label, mat) { console.log(label) console.log(mat.join`\n`) } ;[ {a:[[1,2],[3,4]],b:[[5,6],[7,8]] }, {a:[[1],[2]],b:[[1,2]]}, {a:[[16,2,3,13],[5,11,10,8],[9,7,6,12],[4,14,15,1]],b:[[1,1],[0,1]]}, {a:[[2]],b:[[5]]} ].forEach(t=>{ show('A',t.a) show('B',t.b) show('A⊗B',f(t.a,t.b)) show('B⊗A',f(t.b,t.a)) console.log('-----------------') })

<pre id=O></pre>

Uiua, 10 bytes

/⊂≡/≡⊂פ¤: 

Try it!

I wouldn't be surprised if this was suboptimal. If we can take input reversed, remove : flip for -1 byte. ¤ fixing twice forces × multiply to vectorize, and then we / reduce ⊂ join across various ≡ rows.

JavaScript (Node.js), 49 bytes

f=(a,b)=>a.flatMap(t=>b.map(u=>t.map?f(t,u):t*u)) 

Try it online!

Nothing special. Can save extra byte if [].at exist