# Vediamo in dettagli la struttura dell'algoritmo. Esso è comunque
# incorporato nella libreria accessibile con:
# source("http://macosa.dima.unige.it/r.R")
#
x <- c(1,2,3,5); y <- c(131,113,89,7); n <- length(x)
a <- sum(x); b <- sum(x^2); c <- sum(x^3); d <- sum(x^4)
e <- sum(y); f <- sum(x*y); g <- sum(x*x*y)
ma <- matrix(data = c(n,a,b,a,b,c,b,c,d), nrow = 3, ncol = 3)
noti <- matrix(data = c(e,f,g), nrow = 3, ncol = 1)
S <- solve(ma,noti); S
           [,1]
[1,] 137.327273
[2,]  -1.945455
[3,]  -4.818182
plot(x,y)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
F <- function(x) S[1]+S[2]*x+S[3]*x^2
curve(F, add=TRUE, col="blue")

Prova con questi altri dati:

x <- c(-5,-2.6,-0.4,2.2,3); y <- c(-1.3,-2.6,-1,4,5.4)
Ottieni:
-0.502079825 
1.35599275 
0.236170866 

Se ho solo tre punti, trovo la parabola che passa esattamente per essi:

x <- c(-1,2,3); y <- c(4,-1,6)
Ottengo:
-2 
-3.83333333   (ovvero: -23/6)
13/6