Conceptual basis for formula for sum of n first integers raised to power k, Most programmers (including me) are painfully aware of quadratic behavior resulting from a loop that internally performs 1, 2, 3, 4, 5 and so on operations per iteration,
$$\sum_{i=1}^n i = \frac{n \left(n+1\right)}{2} $$
It’s very easy to derive, e.g. considering the double of the sum like $(1+2+3) + (3+2+1) = 3\times4$ .
For the sum of squares there is a very similar formula,
$$\sum_{i=1}^n i^2 = \frac{n(n+\frac{1}{2})(n+1)}{3}$$
But the only way I found to derive that was to assume it would sum to a cubic polynomial $f(x)=Ax^3+Bx^2+Cx+D$, and solve for $f(n)-f(n-1)=n^2$.
I’m guessing that there is a much simpler system and concept here, generalizing to $\sum_{i=1}^n i^k$?

These "power sum" polynomials are known as Bernoulli polynomials, and have been studied for centuries, and there is a vast literature about them. There are many inductive formulae relating them. As you observed, the right thing to do is to consider a polynomial $f_{k}(x)$ with $f_{k}(0) = 0$
and $f_{k}(x+1) - f_{k}(x) = (x+1)^{k}$, where $k$ is a chosen positive integer.
Notice that is a polnomial of finite degree satisfies this equation for all positive integers, then it must satisfy it for all real $x$, and furthermore, it is unique.
Notice also that, given such a polynomial exists, we must have $f(-1) =0$, so that both $x$ and $x+1$ must be factors
for $f_{k}(x)$ for every $k$, given that $f_{k}$ exists. How do we know that the 
polynomial $f_{k}(x)$ always exists? There are many ways to see this. I like an inductive
approach. The polynomial $f_{1}(x) = \frac{x(x+1)}{2}$ gets us started. How can we find $f_{k+1},$ given $f_{k}$? Well, one way to do it is to notice that if we had $f_{k+1}$
and differentiated its defining equation, we would obtain
$f_{k+1}^{'}(x+1) - f_{k+1}^{'}(x) = (k+1)(x+1)^{k}$, which is nearly the defining 
equation for $f_{k}$, apart from a factor $k+1$ and the possible addition of a constant. Hence, if it is to exist, we should have $f_{k+1}(x) = c(k+1)x + d(k+1) + (k+1)\int_{-1}^{x} f_{k}(t) dt$  for certain constants $c(k+1)$ and $d(k+1)$. We can determine the constants $c(k+1)$ and $d(k+1)$.
Since we need $f_{k+1}(0) = 0$, we must have $d(k+1) = -(k+1)\int_{-1}^{0} f_{k}(t)dt$.
Since we need $f_{k+1}(-1) = 0$, we need $c(k+1) = d(k+1)$.  Hence we have uniquely
specified a polynomial $f_{k+1}$ (of degree $k+2$) with the right properties.
It is $f_{k+1}(x) = -(x+1)(k+1)\int_{-1}^{0}f_{k}(t)dt + (k+1) \int_{-1}^{x} f_{k}(t)dt$.
This can be rewritten as 
$$ \frac{f_{k+1}(x)}{k+1} = x \int_{0}^{-1}f_{k}(t)dt + \int_{0}^{x} f_{k}(t)dt$$
if preferred.
