The following MATLAB code generates the 2D Laplacian matrix using a Kronecker product approach.

function A=A(N)

% Assemble the system matrix A
e = ones(N,1);
D = spdiags([e -2*e e], -1:1, N, N);
I = speye(N);
A = kron(D,I)+kron(I,D);
h = 1/(N+1);
A=A/h^2;

`

I need to convert it to python in order to solve a PDE boundary value problem.

My first attempt involved using SMOP. Ran into a heap of errors and decided it wasn't worth the trouble.

My second attempt involved converting the code directly to python by using the numpy and scipy libraries.

Here's my attempt so far

import numpy as np
import scipy.sparse as sps


def build_laplacian_2D_kron(N):
    e = np.ones(shape=(N, 1))
    data = np.concatenate([e, -2*e, e] , axis= -1)
    diags = np.array([-1, 0, 1])                                                                                                                                                                                                                                                                                     
    D = sps.spdiags(data.transpose(), diags, N, N)
    I = sps.speye(N)
    A = np.kron(D, I) + np.kron(I, D)
    h = 1 / (N + 1)
    A = A / h ** 2
    A = sps.spmatrix.tocsc(A)

    return A


if __name__ == '__main__':
    build_laplacian_2D_kron(6)

When I attempt to run the code the following error shows up

ValueError: offsets array must have rank 1

The error is raised by the function scipy.sparse.spdiags when invoking sps.spdiags on line 9. I presume that the issue is mainly with the way I am initializing data and diags. I have been trying to fix it for quite some time now but to no avail.

Any help would be highly appreciated

Update: With the recommendations from hpaulj I fixed the errors in the python implementation and updated the code above. Here's the output for N = 2

  (0, 0)    -36.0
  (1, 0)    9.0
  (2, 0)    9.0
  (3, 0)    0.0
  (0, 1)    9.0
  (1, 1)    -36.0
  (2, 1)    0.0
  (3, 1)    9.0
  (0, 2)    9.0
  (1, 2)    0.0
  (2, 2)    -36.0
  (3, 2)    9.0
  (0, 3)    0.0
  (1, 3)    9.0
  (2, 3)    9.0
  (3, 3)    -36.0

Process finished with exit code 0

However it does not quite match with the output in matlab(which excludes the 0 entries in the resulting matrix)

ans =

   (1,1)      -36
   (2,1)        9
   (3,1)        9
   (1,2)        9
   (2,2)      -36
   (4,2)        9
   (1,3)        9
   (3,3)      -36
   (4,3)        9
   (2,4)        9
   (3,4)        9
   (4,4)      -36

What modification can I make to my python implementation to exclude the 0 entries in my output? Any help would be appreciated.