using LinearAlgebra, SparseArrays
    N = 10
    D1 = sparse(Tridiagonal(ones(N-1),-2ones(N),ones(N-1)))
    Id = sparse(I,N,N)
    A = kron(D1,Id) + kron(Id,D1)

You just use the Kronecker product to build the matrix in the different directions. This, along with the relationship to convolutional neural networks, is described in the MIT 18.337 course notes: https://mitmath.github.io/18337/lecture14/pdes_and_convoluti...