cubicspline <- function ( x, y )

#*****************************************************************************80
#
## cubicspline computes a cubic spline to data.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    23 February 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
{
  source ( "/home/burkardt/public_html/r_src/tridiagmatrix/tridiagmatrix.R" )

  n <- length ( x )
  avec <- rep ( 0, n - 1 )
  bvec <- rep ( 0, n - 1 )
  dvec <- rep ( 0, n - 1 )
  vec <- rep ( 0, n )
  deltax <- rep ( 0, n - 1 )
  deltay <- rep ( 0, n - 1 )
#
#  Find delta values and the A vector.
#
  for ( i in 1 : ( n - 1 ) )
  {
    avec[i] <- y[i]
    deltax[i] <- x[i+1] - x[i]
    deltay[i] <- y[i+1] - y[i]
  }
#
#  Assemble a tridiagonal matrix of coefficients.
#
  Au <- c ( 0.0, deltax[2:(n-1)] )
  Ad <- c ( 1.0, 2.0 * ( deltax[1:(n-2)] + deltax[2:(n-1)] ), 1.0 )
  Al <- c ( deltax[1:(n-2)], 0.0 )

  vec[0] <- 0
  for ( i in 2 : ( n - 1 ) )
  {
    vec[i] <- 3.0 * ( deltay[i] / deltax[i] - deltay[i-1] / deltax[i-1] )
  }
  vec[n] <- 0

  cvec <- tridiagmatrix ( Al, Ad, Au, vec )
#
#  Compute B and D vectors from the C vector.
#
  for ( i in 1 : ( n - 1 ) )
  {
    bvec[i] <- ( deltay[i] / deltax[i] ) - ( deltax[i] / 3.0 ) * ( 2.0 * cvec[i] + cvec[i+1] )
    dvec[i] <- ( cvec[i+1] - cvec[i] ) / ( 3.0 * deltax[i] )
  }

  return ( list ( a = avec, b = bvec, c = cvec[1:(n-1)], d = dvec ) )
}