Newton’s Divided Difference Polynomial Method of Interpolation

First we define the data we have
x: the values of x
y: the value of the function of x f(x)

x<- c(2,4,5,6,7)
y<- c(3,5,1,6,9)

Then we define the function of the newton interpolation

new_interp <- function(no) {
  n <- length(x)
  mat <- cbind(x, y)
  
  A <- matrix(c(rep(0, ((n) ^ 2 - n))),
              nrow = n,
              ncol = n - 1,
              byrow = T)
  
  for (j in 1:(n - 1))
    for (i in 1:(n - j)) {
      if (j == 1)
        A[i, j] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
      else
        A[i, j] = (A[(i + 1), (j - 1)] - A[(i), (j - 1)]) / (x[i + j] - x[i])
    }
  
  nam <- as.list(seq(1, n - 1, by = 1))
  colnames(A) <- nam
  #print(A)
  
  A1 <- cbind(mat, A)
  Adf <- as.data.frame(A1)
  
  cat(" Divided Difference Table \n ________________________ \n")
  print(Adf)
  
  #no<-2.5
  sum = y[1]
  for (i in 1:(n - 1)) {
    prod = 1
    for (j in 1:i)
      prod = prod * (no - x[j])
    sum = sum + prod * A[1, i]
  }
  #print(sum)
  cat("\n f", n - 1, "(", no, ")= ", sum, "\n")
}

finaly we check the function
Obtain the interpolating polynomials for 2.5

new_interp(2.5)
##  Divided Difference Table 
##  ________________________ 
##   x y  1         2         3      4
## 1 2 3  1 -1.666667  1.541667 -0.675
## 2 4 5 -4  4.500000 -1.833333  0.000
## 3 5 1  5 -1.000000  0.000000  0.000
## 4 6 6  3  0.000000  0.000000  0.000
## 5 7 9  0  0.000000  0.000000  0.000
## 
##  f 4 ( 2.5 )=  12.07031


REFRENCE
Newton’s divided difference formula in R programming language
link: https://www.youtube.com/watch?v=gic5uZBYd-A