The paper gives the rate of convergence you get with Polynomial interpolants in Chebyshev nodes:
> If f has v derivatives, with the vth derivative being of bounded variation V, then ||f - p_n|| = O(V n^{-v}) as n -> ∞
and
> If f is analytic, the convergence is geometric, with ||f - p_n|| = O(p^{-n}) for some p > 1
You will not get that good of a rate of convergence with cubic splines. See https://www.researchgate.net/publication/243095286_On_the_Or...
This is further explained in Trefethen's book https://www.amazon.com/Approximation-Theory-Practice-Applied...
Quoting from Ch 14
> In fact, polynomial interpolants in Chebyshev points are problem-free when evaluated by the barycentric interpolation formula. They have the same behavior as discrete Fourier series for period functions, whose reliability nobody worries about. The introduction of splines is a red herring: the true advantage of splines is not that they converge where polynomials fail to do so, but that they are more easility adapted to irregular point distributions and more localized.
You can see also the software package https://www.chebfun.org/ for Chebyshev interpolations with Matlab and https://github.com/rnburn/bbai for Chebyshev interpolation of arbitrary dimension functions with sparse grids for Python. And here is a quick notebook for an experiment you can run that will compare Chebyshev interpolants to cubic splines: https://github.com/rnburn/bbai/blob/master/example/13-sparse...