This is an excerpt from my math models textbook. It’s about Lagrange Polynomials which is a technique that lets you fit a polynomial to a set of any number of unique points (x_1,y_1) … (x_n,y_n) so long as all your x-values are different (otherwise it wouldn’t be a function, and couldn’t be a polynomial). The polynomial you’ll calculate will be the unique, lowest degree polynomial that passes through all points.
It’s not too bad, once you consider that everything in each term is a constant value, except for “x” itself.
So the numerator of each term is the product of three linear factors, like (x-4)(x-2)(x-6), which should product a cubic, like x³ - 12x² + 44x - 48. Then the denominator of each term is a pure constant, so it would be like taking that cubic and dividing it by 4, getting 0.25x³ - 3x² + 11x - 12. Then the yₙ terms are also constants, so no different than doing something like multiplying by 2, getting you something like 0.5x³ - 6x² + 22x - 24, if I take that example a bit too far. And at that point it’s just the sum of four cubics, which will be cubic as long as the x³ terms don’t perfectly cancel out - which I believe would only happen if the four pairs of points used to make the function had all four points perfectly lying on the same line or parabola.
The construction’s also pretty clever: OP said the point was to fit the function to the four points (x₁, y₁), (x₂, y₂), (x₃, y₃), (x₄, y₄). Let’s say we set x = x₂, then. Because (x-x₂) appears in the numerator of every term but the second term, every term but the second term will have a 0 in its numerator and cancel out - so we only need to consider the second term. Its numerator is then (x₂-x₁)(x₂-x₃)(x₂-x₄) - exactly the same as its denominator. So they both cancel out, leaving only y₂ - meaning we get P₃(x₂) = y₂, as desired.