Spline interpolation is a rich and extensively studied field, dating back to at least the 1940s. Here are a few comments about related work. The interpolation algorithm itself can be adjusted in numerous ways. One effective algorithm for reducing overshooting artifacts is Akima spline interpolation, but there are many others. Cubic splines have noteworthy theoretical properties; you can search for "energy minimizer cubic spline" to find out more. As already mentioned in @creata's comment, instead of using Gaussian elimination, you can take advantage of the fact that the matrix is banded and solve it using a banded solve algorithm (search for DGBSV in LAPACK). There is also research that explores the relationship between shift-invariant spaces and functions that can be expressed as a linear combination of splines. Finally, multivariate interpolation is an extension of the 1-d theory to n-d. This is used in image and surface interpolation, for example.