EGauss<-function(A,b){ # Função para escalonamento de sistemas lineares
n=length(A[,1])        # sem troca de linhas
for (j in 1:(n-1)){
for ( i in (j+1):n){
m= A[i,j]/A[j,j]
A[i,]=A[i,]-m*A[j,]
A[i,j]=0
b[i]=b[i]-m*b[j]
}
}
return(list(R=A,v= b))
}

TSup<-function(R,v){ # Resolução de sistema triangular superior
n=length(v);u=0*v    # Resolve o sistema Ru=v, com R triangular superior
u[n]=v[n]/R[n,n]
for ( i in (n-1):1){
p=v[i]
for (j in (i+1):n){
p=p-u[j]*R[i,j]
}
u[i]=p/R[i,i]
}
u
}

#---------------------------------------------------------------------------
# Problema de interpolação polinomial

n=41
a=-5;b=5;x=0*1:n

for ( k in 1:n){  # Nós de Chebyshev reescalados.
x[k]=(a+b)/2+(b-a)*cos((2*k-1)/(2*n)*pi )/2}

y=1/(1+x^2)

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

n=length(x)
A=matrix(0,n,n); 
for ( j in 1:n){A[,j]=x^(j-1)}


v=y
B=t(A)%*%A; b=t(A)%*%v       

B=EGauss(B,b);R=B$R;v=B$v


u=TSup(R,v); u
print("Teste de erro"); A%*%u-y


poli<-function(s){
    p=u[1]
    for ( i in 2:n){p=p+u[i]*(s)^(i-1)}
    return(p)
} 

curve(poli,-5,5,col="blue",ylim=c(0,1))
points(x,y,col="red")        # Teste de ajuste.


Erro<-function(s){1/(1+s^2)-poli(s)}
curve(Erro,-5,5)