Lagrange polynomial interpolation

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

x = np.array([-2, -1, 0, 1, 2, 3])
y = np.array([-4, 4, -2, -3, 3, -5])
plt.plot(x, y);

M = len(x)
w = y

fig, axs = plt.subplots(M, M)

p = sp.poly1d(0.0)
for j in range(M):
    pt = sp.poly1d(w[j])
    for k in range(M):
        ax = axs[j][k]
        ax.axis('off')
        if k == j:
            ax.plot(x, y, 'r')
            continue
        fac = x[j]-x[k]
        pt *= sp.poly1d([1.0, -x[k]])/fac
        ax.plot(pt(x), 'k')

    xp = np.linspace(x.min(), x.max(), 100)
    axs[j][j].plot(xp, p(xp), 'g')
    p += pt

plt.tight_layout()