Polynomial interpolation in python

Question

I am studying function approximation and while trying to understand/implement polynomial interpolation I've found an example here. I find the code below a good example to understand what is actually going on instead of using ready functions, however it doesn't run:

Defining the interpolation algorithm. Essentially, we are trying to come up with a representation of true f as a linear combination of basis functions (psi-s).

import sympy as sym
def interpolation(f, psi, points):
    N = len(psi) - 1 #order of the approximant polynomial
    A = sym.zeros((N+1, N+1)) # initiating the square matrix, whose regular element is psi evaluated at each nodes
    b = sym.zeros((N+1, 1)) # original function f evaluated at the selected nodes
    psi_sym = psi  # save symbolic expression
    # Turn psi and f into Python functions
    x = sym.Symbol('x')
    psi = [sym.lambdify([x], psi[i]) for i in range(N+1)]
    f = sym.lambdify([x], f)
    for i in range(N+1):
        for j in range(N+1):
            A[i,j] = psi[j](points[i])
        b[i,0] = f(points[i])
    c = A.LUsolve(b) #finding the accurate weights for each psi
    # c is a sympy Matrix object, turn to list
    c = [sym.simplify(c[i,0]) for i in range(c.shape[0])]
    u = sym.simplify(sum(c[i,0]*psi_sym[i] for i in range(N+1)))
    return u, c 

True function f:= 10(x-1)^2 -1, nodes: x0:= 1 + 1/3 and x1:= 1 + 2/3. Interval: [1,2].

x = sym.Symbol('x')
f = 10*(x-1)**2 - 1
psi = [1, x] # approximant polynomial of order 1 (linear approximation)
Omega = [1, 2] #interval
points = [1 + sym.Rational(1,3), 1 + sym.Rational(2,3)]
u, c = interpolation(f, psi, points)
comparison_plot(f, u, Omega) 

The code doesn't run. The error occurs in line

A = sym.zeros((N+1, N+1))  

Error message: ValueError: (2, 2) is not an integer
But A isn't supposed to be an integer, it is a square matrix whose each element is psi evaluated at each node. f = A*c.