How to do polynomial interpolation given a set of points?

Question

Given x points, I know that the interpolated polynomial is at most of degree x-1. However, I do not know how to generate such a polynomial (not a piecewise interpolation). I tried via the following code, but this does not seem to work:

def interpolate(xis, fis, phis):
    A = np.asmatrix([[phi(xi) for phi in phis] for xi in xis])
    weights = np.linalg.solve(A, fis)
    def y(x):
        return sum(map(lambda w, phi: w*phi(x), weights, phis))
    return y

def bf(x_i):
    def basis(x_j):  
        if x_i == x_j: 
            return 1 
        else: 
            return 0
    return basis 
bs = list(map(bf, xis))
f = interpolate(xis, fis, bs)
plot_function_and_data_points(f,xis,fis)

plot:

Does anyone see what I'm doing wrong here?