x=np.array([0,1,2,4,5,6,7,9,10,11])
y=np.array([0,0,1,1,0,0,1,1,0,0])
n=len(x)
h=np.zeros(n-1)
for i in range (0,n-1):
    h[i]=x[i+1]-x[i]

A=np.zeros([n,n])

for i in range (n-2):
    A[i,i]=h[i]
    A[i,i+1]=2*(h[i]+h[i+1])
    A[i,i+2]=h[i+1]
A[-2,0]=1
A[-1,-1]=1


b=np.zeros([n])
for i in range (n-2):
    b[i]=6*(((y[i+2]-y[i+1])/(h[i+1]))-((y[i+1]-y[i])/h[i]))
b[-2]=0
b[-1]=0

M=np.linalg.solve(A,b)


a=np.zeros(n-1)
b=np.zeros(n-1)
c=np.zeros(n-1)
d=np.zeros(n-1)
for i in range (n-1):
    a[i]=(M[i+1]-M[i])/(6*h[i])
    b[i]=M[i]/2
    c[i]=((y[i+1]-y[i])/h[i])-((h[i]/6)*(M[i+1]+2*M[i]))
    d[i]=y[i]



def interpol (xin,x,y):
    n=len(x)
    elem=0
    for i in range (n):
        g=1
        for j in range (n):
            if i!=j:
                g*= (xin-x[j])/(x[i]-x[j])
        elem+= g*y[i]
    return elem
        