def n(j,axis,x):
n=1
for i in np.arange(j):
n *=(axis-x[i])
return n
def divdif(y,x,i,n):
if n==0:
return y[i]
if n==1:
return (y[i]-y[i-1])/(x[i-1]-x[i])
else:
return (divdif(y,x,i,n-1)-divdif(y,x,i-1,n-1))/(x[i-n]-x[i])
def N(axis,x,y):
N=0.
for j in range(len(x)):
N +=divdif(f(x),x,-1,j)*n(j,axis,x)
return N
for order in dano:
x=np.linspace(-1.,1.,order,endpoint=True)
plt.plot(x,f(x),'d')
plt.plot(xaxis,N(xaxis,x,f(x)),label='n='+str(order))