A tiny 1D Laplacian is sparse. spsolve solves Ax = b without a dense inverse.

Goal

Solve a 5×5 tridiagonal system and print x.

from scipy import sparse
from scipy.sparse.linalg import spsolve
n = 5
diags = [np.ones(n - 1), -2 * np.ones(n), np.ones(n - 1)]
A = sparse.diags(diags, [-1, 0, 1], format='csr')
b = np.ones(n)
x = spsolve(A, b)
print(x.round(4))
print('residual', np.max(np.abs(A @ x - b)))

from scipy import sparse
A = sparse.eye(4, format='csr') * 2
print(A.toarray())

from scipy.sparse.linalg import spsolve
from scipy import sparse
A = sparse.diags([1, 1, 1], format='csr')
print(spsolve(A, np.array([4.0, 5, 6])))

from scipy import sparse
print(sparse.kron(sparse.eye(2), sparse.eye(2)).toarray())