Kronecker product and finite difference discretization for poisson equation

Question

In this notebook from MIT's Intro to Linear PDEs course, it is unclear to me why the Kronecker product is used to formulate the coefficients matrix $A$ for solving the linear system of equations $A u = b$. Given Poisson's equation on the unit square with homogeneous boundary conditions (see Finite Difference Methods for Poisson Equation by Long Chen, UCI Math for general treatment of Neumann and Dirichlet boundary conditions), the goal is to find $u(x,y)$ satisfying

$$ \nabla^2u = f. $$

The second-order centered difference approximation is represented by a matrix $D \in \mathbb{R}^{n \times n}$ for an $n \times n$ equispaced computational grid, but why is the Kronecker product and identity matrix used in $A = I\ \otimes D + D\ \otimes I$ to form the coefficient matrix? Clearly, the ordering of the arguments for the Kronecker product also determines which dimension ($x$ or $y$) the second-order centered difference approximation operates on (from the MIT notes, the matrix form of the equation is $UD + DU = F$ and the matrix $D$ is a linear operator on $U$), but what about for the 3D case?

Here is the Julia code from the notebook that builds the coefficient matrix $A$ for reference:

using LinearAlgebra
using SparseArrays

## build n x n 2nd central difference matrix
function cdiff2(n)
    h = 2 / (n+1)
    D = spdiagm(-1 => ones(n-1), 0 => -2*ones(n), 1 => ones(n-1))
    D = D / h^2
end

## build n^2 x n^2 discretization of Laplace operator with zero Dirichlet BCs
function lap(n)

    # sparse second difference matrix and identity
    D = cdiff2(n)
    spI = sparse(I,n,n)

    # Kronecker product
    A = kron(spI,D)+kron(D,spI)
end

Edit:

A related question/answer is here