You don't need to use O(n^3) Gaussian elimination to find natural cubic splines: you can write the system of equations as a tridiagonal matrix, which can be solved in linear time.
You don’t need to solve a system of equations at all, and it’s unusual to do so unless you have very specific requirements, right? All of the reasons the article used to justify going that direction can be met (more or less, depending on more nuance than the article used) with existing spline types. Plus you can trade some of the overshoot properties of the article’s solution for just a little bit of extra curvature, with a Catmull-Rom spline for example, which might be preferable for many people for multiple reasons.
https://splines.readthedocs.io/en/latest/euclidean/natural-u...