def solve(n):
    nrange = np.arange(n)
diag = np.ones(n - 1)
n_to_M = spdiags(-2.0 * diag, 0, n - 1, n - 1)
n_to_N = spdiags([n * diag, -nrange[-1:0:-1]], [0, 1], n - 1, n - 1)
m_to_M = lil_matrix(n_to_N)
m_to_M[1:, (0)] = -nrange[1:-1].reshape((n - 2, 1))
m_to_N = lil_matrix((n - 1, n - 1))
m_to_N[:, (0)] = -nrange[1:].reshape((n - 1, 1))
coeff_mat = hstack((vstack((n_to_M, n_to_N)), vstack((m_to_M, m_to_N))))
const = n * np.ones((2 * (n - 1), 1))
return spsolve(coeff_mat.tocsr(), const).reshape((-1, 1))
Using scipy sparse matrices to solve system of equations