x=c(9,8,1,6,10,5,3);x; y=c(3,2,1,5,12,7,13);y

plot(x,y,col="red")

n=length(x)

# Diferenças divididas

A=matrix(0,n,n);
A[,1]=y

for (j in 2:n){
for ( i in 1:(n-j+1)){
A[i,j]=(A[i,j-1]-A[i+1,j-1])/(x[i]-x[i+j-1])
}
}

A

N<-function(i,s){
p=1
for (j in 1:i){p=p*(s-x[j])}
p
}

poli<-function(s){ 
q=A[1,1]
for (i in 2:n){q=q+A[1,i]*N(i-1,s)}
q}

curve(poli,min(x),max(x))
points(x,y,col="red")