# splint.r          cubic spline interpolation
 #
  npts = 100                        # points for plotting
   n = 7                            # how many points for interpolation
   xi <- 4*(2*(1:n)-n-1)/(n-1)      # equally spaced interpolation pts    
   z <- 1/(1+xi*xi)                 # function values
   szx <-splinefun(xi,z,method="natural")# interp cubic spline function
 # now get ready to plot
   x <- 4*(2*(1:npts)-npts-1)/(npts-1)
   zx <- 1/(1+x*x)                   # true function
   s <- szx(x,deriv=0)               # interpolated values
   sp <- szx(x,deriv=1)              # deriv of spline
   cbind(x,zx,s,sp)
   plot(c(x,x),c(zx,s))              # Figure 7.3
 # clean up & go
  rm(list=ls())
  q()