import numpy as np
from scipy.sparse import diags_array, issparse
from scipy.sparse.linalg import spsolve

n = 5
A = diags_array(
    [-np.ones(n - 1), 4 * np.ones(n), -np.ones(n - 1)],
    offsets=[-1, 0, 1],
    shape=(n, n),
    format="csr",
)
b = np.array([2.0, 1.0, 0.0, 1.0, 2.0])

x = spsolve(A, b)
residual = A @ x - b

print(f"matrix format: {A.format}")
print(f"matrix shape: {A.shape}")
print(f"stored values: {A.nnz}")
print("solution:", np.round(x, 6))
print(f"residual norm: {np.linalg.norm(residual):.3e}")
print("passes tolerance:", np.linalg.norm(residual) < 1e-12)
print("matrix is sparse:", issparse(A))