x=0*seq(-5,5,by=0.25); n=length(x)

a=-5;b=5
for ( k in 1:n){
x[k]=(a+b)/2+(b-a)*cos((2*k-1)/(2*n)*pi )/2}
y=1/(1+x^2)

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])
}
}

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")