It seems more can be said about the system,
$x+y+z =1$
$x^3 +y^3 +z^3 +xyz = x^4 + y^4 +z^4 +1$
which, after an adjustment, can be connected to $p=7$ and the heptagon. It is just 2 equations in 3 unknowns, hence is under-determined and should have infinitely many solutions. Interestingly, it seems there are two classes of solutions: 1st, using only square roots of square roots (quartics) and 2nd, roots of cubics.
I. Roots of quartics
Solve for $x$ in the first equation, then $y$ in the second, and we have the solution,
$$x = \frac{1-z+\sqrt{-3z^2+2z+\sqrt{-7}}}2\\
y = \frac{1-z-\sqrt{-3z^2+2z+\sqrt{-7}}}2$$
where we can choose either sign of $\pm\sqrt{-7}$, with free parameter $z$. Thus, the system has infinitely many solutions where one variable can be real ($z$), but the other two are complex.
II. Roots of cubics
As pointed out by Thomas Andrews, we can use the elementary symmetric polynomials. Rephrase the system as,
$x+y+z =1 $
$x^3 +y^3 +z^3 +xyz = k$
$x^4 + y^4 +z^4 +1 = k$
and we now have 3 equations in 3 unknowns for some free parameter $k$. The $(x,y,z)$ then are the three roots of the cubic,
$$16 x^3 - 16x^2 + 4x - 4k + 1 = (4 x - 3)\sqrt{-7}$$
Because of the $\sqrt{-7},\,$ this generally yields only complex roots. For example, let $k=2$,
$$16 x^3 - 16x^2 + 4x - 7 = (4 x - 3)\sqrt{-7}$$
and the cubic has roots,
$$(x,y,z) \approx (-0.4134 - 0.7567 i,\; 0.3171+ 0.6174 i,\; 1.0962 + 0.1393 i)$$
which solve the OP's system. But they are all complex in contrast to the previous method where one of $(x,y,z)$ can be real.
III. Real roots
However, there are similar systems with all $(x,y,z)$ solvable in the reals (as the OP wished) and where the infinite number of solutions are again of two kinds: constructible numbers (which uses only square roots), and non-constructible numbers (which uses higher roots). Consider the system,
$x+y+z =1$
$x^3 +y^3 +z^3 +xyz = x^4 + y^4 +z^4 \color{red}{-10}$
which differs from the OP's only in the red number. Using the same approach as above, then:
A. First class (constructible)
The $(x,y,z)$ are,
$$x = \frac{1-z+\sqrt{-3z^2+2z+9}}2\\
y = \frac{1-z-\sqrt{-3z^2+2z+9}}2$$
and we can choose infinitely many rational $z$ such that $-3z^2+2z+9\geq0$.
B. Second class (non-constructible)
The $(x,y,z)$ are the three roots of the cubic,
$$x^3 - x^2 - 2x - n + 1 = 0$$
and we can choose infinitely many real $n$ such that $-27n^2 + 14n + 49\geq0$ so all roots are real. The most well-known example is $n=0$ and the cubic,
$$x^3 - x^2 - 2x + 1 = 0$$
whose roots are,
$$(x,y,z)= \big({-2}\cos(2\pi/7),\;-2\cos(4\pi/7),\;-2\cos(6\pi/7)\big)$$
which appear in the heptagon and solve the given system, namely,
$x+y+z =1$
$x^3 +y^3 +z^3 +xyz = x^4 + y^4 +z^4 \color{red}{-10}$
And so on for other $n$. Thus, if this system had been the example in Arthur Engel's problem solving strategies, then both solution classes would have to be mentioned.