Formula for $1^k+2^k+3^k...n^k$ for $n,k \in \mathbb{N}$

https://en.m.wikipedia.org/wiki/Faulhaber%27s_formula seems to deliver exactly what i needed. Original answer by MooS

I'll give you a derivation which I don't think is all that known. You just need to know four things:

(1) A version of the umbral Taylor series

Suppose that we have a sequence:

$$a_0,a_1,...a_k$$

And we want to find the function of $n$ that defines $a_n$.

To do this we start by letting $a_{n+1}-a_n=\Delta a_n$ and we call this operation on $a_n$ the forward difference. Then given $\Delta a_n$ we can find $a_n$. Sum both sides of the equation from $n=0$ to $x-1$, and note that we have a telescoping series:

$$\sum_{n=0}^{x-1} \Delta a_n=\sum_{n=0}^{x-1} (a_{n+1}-a_n)=a_{x}-a_{0}$$

Hence $a_n=a_0+\sum_{i=0}^{n-1} \Delta a_i$. Also $\Delta a_n=\Delta (0)+\sum_{i=0}^{n-1} \Delta^2 a_n$...and so forth. Using this we must have if the series converges:

$$a_n=a_0+\Delta (0) \sum_{x_0=0}^{n-1} 1+\Delta \Delta (0) \sum_{x_0=0}^{n-1} \sum_{x_1=0}^{x_0-1} 1+\Delta \Delta \Delta (0) \sum_{x_0=0}^{n-1} \sum_{x_1=0}^{x_0-1} \sum_{x_2=0}^{x_1-1} 1+\cdots$$

Where $\Delta^i (0)$ denotes the first term ($n=0$) of the $i$ th difference sequence of $a_n$.

Through a combinational argument, If we take $\Delta^0 (0)=a_0$ and ${n \choose 0}=1$ we may get:

$$a_n=\sum_{i=0}^{\infty} \Delta^i(0) {n \choose i}$$

(2) Hockey stick identity Proof of the Hockey-Stick Identity: $\sum\limits_{t=0}^n \binom tk = \binom{n+1}{k+1}$

(3) The forward difference of a polynomial of degree $k$ is of degree $k-1$ (follows from binomial theorem) hence $k+1$ forward differences survey results in $0$ as $k$ forward differences must result in a constant sequence.

(4) A polynomial of degree $k$ can be uniquely defined by the sequence $f(0),f(1),f(2),...f(k)$.

Using (1), (2), (3), and (4) you may come up with the formula:

$$\sum_{s=1}^{n} s^k=\sum_{s=1}^{k} b_s{n+1\choose s+1}$$

By considering how to represent $a_s=s^k$ as a sum of multiples of binomials.

$$s^k=0+b_1{ s \choose 1}+b_2{s \choose 2}+\cdots+ b_k{s \choose k}+0+0+\cdots=\sum_{i=1}^{k} b_i{s \choose i}$$.
$$\sum_{s=1}^{n}\sum_{i=1}^{k} b_i{s \choose i}=\sum_{s=1}^{n}b_1{ s \choose 1}+\sum_{s=1}^{n}b_2{s \choose 2}+\cdots +\sum_{s=1}^{n}b_k{s \choose k}$$.
Taking the constants $b_1,b_2,...$ out of the sums and utilizing the hockey stick identity we get our desired result.

Now see Function $f: \mathbb{Z} \to \mathbb{Z}^n$ related to $\sum_{k=1}^{x} k^n$.

We have that:

$$b_s=s!S(k,s)$$

Where ! denotes the factorial and $S( , )$ denotes the stirling numbers of the second kind.

Hence we have that:

$$\sum_{s=1}^{n} s^k=\sum_{s=1}^{k} s!S(k,s){n+1 \choose s+1}$$

Here are some examples.

$k=1$

$$\sum_{s=1}^{n} s^1=1{n+1 \choose 2}$$

$$=1\frac{(n+1)(n)}{2!}$$

$k=2$

$$\sum_{s=1}^{n} s^2=1{n+1 \choose 2}+2{n+1 \choose 3}$$

$$=1\frac{(n+1)(n)}{2!}+2\frac{(n+1)(n)(n-1)}{3!}$$

$k=3$:

$$\sum_{s=1}^{n} s^3=1{n+1 \choose 2}+6{n+1 \choose 3}+6{n+1 \choose 4}$$

$$=1\frac{(n+1)(n)}{2!}+6\frac{(n+1)(n)(n-1)}{3!}+6\frac{(n+1)(n)(n-1)(n-2)}{4!}$$

Another way to derive a formula for $$S(k,n)=\sum_{s=1}^{n}{s^k}$$ is to use the binomial expansion. To do this, you can start by changing the summation index $s$ to $t+1$. That is:$$S(k,n)=\sum_{t=0}^{n-1}{(t+1)^k}=1+\sum_{t=1}^{n-1}{(t+1)^k}$$Using binomial expansion we have:$$S(k,n)=1+\sum_{t=1}^{n-1}{\sum_{j=0}^{k}{\binom{k}{j}t^j}}$$Changing the order of summation (notice the independence) gets:$$S(k,n)=1+\sum_{j=0}^{k}{\sum_{t=1}^{n-1}{\binom{k}{j}t^j}}=1+\sum_{j=0}^{k}{\left(\binom{k}{j}\sum_{t=1}^{n-1}{t^j}\right)}$$The inner sum equals $S(j,n-1)$, so: $$S(k,n)=1+\sum_{j=0}^{k}{\binom{k}{j}S(j,n-1)}$$Now, by excluding the last two terms of the summation (the terms obtained by $j=k-1$ and $j=k$) we finally get: $$S(k,n)=1+S(k,n-1)+kS(k-1,n-1)+\sum_{j=0}^{k-2}{\binom{k}{j}S(j,n-1)}$$ Since $S(k,n)-S(k,n-1)=n^k$: $$S(k-1,n-1)=\frac{1}{k}\left(n^k-1-\sum_{j=0}^{k-2}{\binom{k}{j}S(j,n-1)}\right)$$ The formula shows a recursive relation between $S(k-1,n-1)$ and lower sums, i.e., $S(0,n-1),S(1,n-1),\dots,S(k-2,n-1)$. By changing $k$ to $k+1$ and $n$ to $n+1$ the original format emerges: $$S(k,n)=\frac{1}{k+1}\left((n+1)^{k+1}-1-\sum_{j=0}^{k-1}{\binom{k+1}{j}S(j,n)}\right)$$ You can now start from $S(0,n)=n$ and derive $S(k,n)$ for $k=1,2,\dots$. For example: $$S(1,n)=\frac{1}{2}\left( (n+1)^2-1-\binom{2}{0}S(0,n)\right)=\frac{1}{2}\left( n^2+2n+1-1-n \right)=\frac{1}{2}n(n+1)$$