R function for performing Neville's Algorithm for approximating an interpolated polynomial at a given point

nevilles_algorithm.R

# Performs Neville's Algorithm (https://en.wikipedia.org/wiki/Neville%27s_algorithm) for 
# interpolating a polynomial given a set of x and y values.
# Parameters:
#   x: x values of interpolating points
#   y: values corresponding to x values
#   x0: desired point to approximate interpolated polynomial
# Returns:
#   List containing approximated value of interpolated polynomial and matrix representing
#   the intermediate values of Q


poly.neville <- function(x, y, x0) {
  
  if (length(x) != length(y)) {
    stop('x and y must be equal length')
  }
  
  n <- length(x)
  q <- matrix(data = 0, n, n)
  q[,1] <- y
  
  for (i in 2:n) {
    for (j in i:n) {
      q[j,i] <- ((x0 - x[j-i+1]) * q[j,i-1] - (x0 - x[j]) * q[j-1,i-1]) / (x[j] - x[j-i+1])
    }
  }
  
  res <- list('Approximated value'=q[n,n], 'Neville iterations table'=q)
  return(res)
}