%matplotlib inline
import pylab as pl
import numpy as np

M, N = 20, 40
np.random.seed(42)
m = np.random.rand(M, N)
m[1:-1, 1:-1] = 0

count = 0
while True:
    count += 1
    avg = 0.25 * (m[:-2, 1:-1] + m[2:, 1:-1] + m[1:-1, :-2] + m[1:-1, 2:])
    err = avg - m[1:-1, 1:-1]
    if np.max(np.abs(err)) < 1e-15:
        break
    m[1:-1, 1:-1] = avg
print count

pl.imshow(m, interpolation="nearest");

from scipy import sparse
from scipy.sparse import linalg

idx = np.arange(M * N).reshape(M, N)

A = sparse.dok_matrix((M*N, M*N))
b = sparse.dok_matrix((M*N, 1))
i = 0
for row in range(M):
    for col in range(N):
        if row == 0 or col == 0 or row == M-1 or col == N-1:
            A[i, idx[row, col]] = 1
            b[i, 0] = m[row, col]
        else:
            A[i, idx[row, col]] = 4
            A[i, idx[row-1, col]] = -1
            A[i, idx[row+1, col]] = -1
            A[i, idx[row, col-1]] = -1
            A[i, idx[row, col+1]] = -1
        i += 1
        

x = linalg.spsolve(A, b)
x.shape = M, N

np.allclose(x, m)