import numpy as np 
import matplotlib.pyplot as plt
def dif_divise(x,y):
    n=len(y)
    coeff=np.zeros([n,n])
    coeff[:,0]=y
    for j in range(1,n):
        for i in range(n-j):
            coeff[i,j]=(coeff[i+1,j-1]-coeff[i,j-1])/(x[i+j]-x[i])
    return coeff
def polynome(x,y,t):
    d=dif_divise(x,y)[0,:]
    n=len(d)-1
    p=d[n]
    for k in range(1,n+1):
        p=d[n-k]+(t-x[n-k])*p
    return p

def f(x):
    return np.cos(x)
x=np.linspace(-5,5,6)
y=[f(t) for t in x]
plt.plot(x,y,'r*')
x1=np.linspace(-5,5,100)
y1=[f(t) for t in x1]
plt.plot(x1,y1,'r')
x0=np.linspace(-5,5,300)
y0=[polynome(x, y, h) for h in x0]
plt.plot(x0,y0,'b')
plt.show()