def pts_interp_reg(f,n):
    L = []
    for j in range(n+1):
       x = -1 + 2*j/n
       L.append( (x,f(x)) )
    return L

RRx.<x> = RR[]
# ou encore
RRx.<x>=PolynomialRing(RealField())

def interp_reg(f,n):
    L = pts_interp_reg(f,n)
    p = RRx.lagrange_polynomial(L)
    return p

f = lambda x: 1/(1+25*x^2)

n=5;plot(interp_reg(f,n)-f(x),-1,1)

a=animate([
  plot([interp_reg(f,n)-f(x),n],-1,1)
  for n in range(5,25)],xmin=-1,xmax=1,ymin=-25,ymax=25)
a.show(delay=100)