There is a question I have based on the fact: If you take a quadratic polynomial with integer coefficients, take the set $\{1,2,3,4,5,6,7,8\}$, make a partition $A=\{1,4,6,7\}$, $B=\{2,3,5,8\}$, and then evaluate the polynomial with the elements of each set you got an equality.
If now with a cubic polynomial and the set $\{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16\}$, how to pick a suitable partition that accomplish the same fact? I mean is there a general formula for bigger degree?
If I understood the question correctly, the following construction will work. We begin with two sets $A_1=\{1,4\}$ and $B_1=\{2,3\}$ (we might equally well similarly partition the set $\{0,1,2,3\}$, but you seem to asking about this set). We observe that $$ \sum_{x\in A_1}x=1+4=5=2+3=\sum_{x\in B_1}x. $$ Because both sets have two elements (or, equivalently, $\sum_{x\in A_1}x^0=\sum_{x\in B_1}x^0$) it follows that for any linear polynomial $f$ we also have $$ \sum_{x\in A_1}f(x)=\sum_{x\in B_1}f(x). $$ We want a similar result to hold for quadratic polynomials $f$, too. To that end we must use larger sets and define $$ A_2=A_1\cup\{4+x\mid x\in B_1\},\qquad B_2=B_1\cup\{4+x\mid x\in A_1\}. $$ We notice that $A_2\cap B_2=\emptyset$ and $A_2\cup B_2=\{1,2,\ldots,8\}.$ We then get the property described in the question: $$ \sum_{x\in A_2}x^k=\sum_{x\in B_2}x^k $$ for all $k=0,1,2$. This is more or less immediate for $k=0,1$. When $k=2$ it follows from the calculation $$ \begin{aligned} \sum_{x\in A_2}x^2&=\sum_{x\in A_1}x^2+\sum_{x\in B_1}(4+x)^2\\ &=\sum_{x\in A_1\cup B_1}x^2+\sum_{x\in B_1}(16+8x)\\ &=\sum_{x\in A_1\cup B_1}x^2+16|B_1|+8\sum_{x\in B_1}x, \end{aligned} $$ and the corresponding calculation $$ \begin{aligned} \sum_{x\in B_2}x^2&=\sum_{x\in B_1}x^2+\sum_{x\in A_1}(4+x)^2\\ &=\sum_{x\in A_1\cup B_1}x^2+\sum_{x\in A_1}(16+8x)\\ &=\sum_{x\in A_1\cup B_1}x^2+16|A_1|+8\sum_{x\in A_1}x. \end{aligned} $$ From our earlier results about sums over sets $A_1$ and $B_1$ it follows that these sums are equal. It then follows that the sums of values of any quadratic polynomials over the sets $A_2$ and $B_2$ agree.
You knew all this, so let's move on. We shall recursively defined two sets $A_k$ and $B_k$ such that $A_k\cap B_k=\emptyset$, $A_k\cup B_k=\{1,2,\ldots,2^{k+1}\}$, and that the power sums are equal $$ \sum_{x\in A_k}x^j=\sum_{x\in B_k}x^j $$ for all $j=0,1,2,\ldots, k$.
Assume that we have succeeded in doing the above up to som value of $k$. We then define $$ A_{k+1}=A_k\cup\{2^{k+1}+x\mid x\in B_k\},\qquad B_{k+1}=B_k\cup\{2^{k+1}+x\mid x\in A_k\}. $$ I claim that this works. Notice that $A_{k+1}$ contains all of $A_k$ and then some new integers in the range $[2^{k+1}+1,2^{k+2}]$. Similarly for $B_{k+1}$. The integers in this new range are equally divided in the two parts, because by the induction hypothesis they were so divided in the previous step. The power sum equality is proven in the same way. We need the binomial formula. It will introduce lower degree power sums, but these were all covered by the induction hypothesis, so we will be fine. Let's roll, and pick an exponent $j$ in the range $0\le j\le {k+1}$. $$ \begin{aligned} \sum_{x\in A_{k+1}}x^j&=\sum_{x\in A_k}x^j+\sum_{x\in B_k}(2^{k+1}+x)^j\\ &=\sum_{x\in A_k}x^j+\sum_{x\in B_k}\left(\sum_{\ell=0}^j{j\choose\ell}2^{(k+1)(j-\ell)}x^\ell\right)\\ &=\sum_{x\in A_k\cup B_k}x^j+\sum_{x\in B_k}\left(\sum_{\ell=0}^{j-1}{j\choose\ell}2^{(k+1)(j-\ell)}x^\ell\right)\\ \end{aligned} $$ Here in the last line $\ell<j=k+1$, so by the induction hypothesis we can replace $B_k$ with $A_k$, and continue $$ \begin{aligned} &=\sum_{x\in A_k\cup B_k}x^j+\sum_{x\in A_k}\left(\sum_{\ell=0}^{j-1}{j\choose\ell}2^{(k+1)(j-\ell)}x^\ell\right)\\ &=\sum_{x\in B_k}x^j+\sum_{x\in A_k}\left(\sum_{\ell=0}^{j}{j\choose\ell}2^{(k+1)(j-\ell)}x^\ell\right)\\ &=\sum_{x\in B_k}x^j+\sum_{x\in A_k}(2^{k+1}+x)^j\\ &=\sum_{x\in B_{k+1}}x^j, \end{aligned} $$ as claimed. The inductive step is thus valid, and claim proven for all $k$.
You were particularly interested in the sets $A_3$ and $B_3$. Here they are $$ A_3=\{1,4,6,7,10,11,13,16\},\qquad B_3=\{2,3,5,8,9,12,14,15\}. $$ The set $A_3$ is the union of $A_2$ and all the numbers gotten by adding $8$ to the elements of $B_2$. Similarly $B_3$ contains all the numbers of $B_2$ as well as those gotten by adding $8$ to an element of $A_2$.
In all the steps knowing that the monomial sums all agree up to a certain degree, of course, implies that the sums of values of any polynomial up to that degree then also agree.
Here is another way to look at the problem. Assume the notation,
$$[a_1,a_2,\dots,a_m]^k = a_1^k + a_2^k +\dots +a_m^k$$
Theorem: (Tarry-Escott) If,
$$[a_1,a_2,\dots,a_m]^k = [b_1,b_2,\dots,b_m]^k,\;\; \text{for}\, k = 1,2,3,\dots,n$$
then for any constant $c$,
$$[a_1,\dots,a_m,\,c+b_1,\dots, c+b_m]^k = [b_1,\dots,b_m,\,c+a_1,\dots, c+a_m]^k \\ \text{for}\, k=1,2,3,\dots, n+1$$
This doubles the number of terms, but you can climb the ladder of powers as high as you like. For example, starting with,
$$[1,4]^k = [2,3]^k,\;\; \text{for}\, k = 1$$
using the theorem, we get,
$$[1,\,4,\,c_1+2,\,c_1+3]^k = [2,\,3,\,c_1+1,\,c_1+4]^k,\\ \text{for}\, k = 1,2$$
true for any $c_1$. If we use $c_1 = 4$, we recover your $[1,\,4,\,6,\,7]^k = [2,\,3,\,5,\,8]^k$. Using the theorem again, we have,
$$[1,\,4,\,6,\,7,\,c_2+2,\,c_2+3,\,c_2+5,\,c_2+8]^k = [2,\,3,\,5,\,8,\,c_2+1,\,c_2+4,\,c_2+6,\,c_2+7]^k,\\ \text{for}\, k = 1,2,3$$
true for any $c_2$, and where we can choose $c_2 = 8$ so terms are from $1,2,\dots, 16.$ And so on for any $k = 1,2,3,\dots,n$.
In general, the statement below holds. A proof of this result uses the corollary in my post here. Take $x=n$ and $x_i=m^i$ for $i=0,1,2,\ldots,r$, and then utilize the fact that the every polynomial $f(T)\in\mathbb{C}[T]$ vanishing at $T=\omega_m^q$ for all $q=1,2,\ldots,m-1$ is divisible by the polynomial $1+T+T^2+\ldots+T^{m-1}$.
Let $m$, $n$, and $r$ be integers such that $m>1$ and $r\geq 0$. The set $$\left\{n,n+1,n+2,\ldots,n+m^r-1\right\}$$ can be partitioned into $m$ subsets $X_1,X_2,\ldots,X_m$ such that $$\sum_{i\in X_1}\,i^l=\sum_{i\in X_2}\,i^l=\ldots=\sum_{i\in X_m}\,i^l$$ for all $l=0,1,2,\ldots,r$. Here, $0^0$ is defined to be $1$.
This can be done by taking $X_j$ to be the set of $n+\lambda$, where the sum of the $m$-ary digits of $\lambda\in\left\{0,1,2,\ldots,p^r-1\right\}$ (in its $m$-ary representation) is congruent to $j$ modulo $m$. It is an interesting question whether the sets $X_j$'s satisfying this property are unique (up to reordering).
