x=seq(-5,5,by=1);x; y=1/(1+x^2);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")

f<-function(s){1/(1+s^2)}

curve(f,-5,5)
points(x,y,col="red")

Erro<-function(s){f(s)-poli(s)}

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