Here is a natural way to solve the system. Continue with,
$$x + y = (x + y)(x^2 − xy + y^2),\>\>\>\>\>
x − y = 2(x − y)(x^2 + xy + y^2)$$
or,
$$(x + y)(x^2 − xy + y^2-1)=(x − y)(x^2 + xy + y^2-\frac12)=0$$
So, there are four cases to be considered,
Case 1: Substitute $x-y=0$ into one of the original equations, say, $y + 3x = 4x^3$, which leads to,
$$4x(x^2-1)=0\implies x= 0, \pm 1$$
Case 2: Similarly, $x+y=0$ leads to,
$$2x(2x^2-1)=0\implies x = 0,\pm \frac 1{\sqrt2}$$
Case 3: Substitute $x^2 − xy + y^2-1=0$ into $y+3x=4x^3$, leading to
$$16x^6-28x^4+13x^2-1=0,\>\> \text{or}$$
$$(x^2-1)(16x^4-12x^2+1)=0\implies x=\pm 1,\pm\frac{\sqrt5\pm 1}{4}$$
Case 4: Similarly $x^2 + xy + y^2-\frac12=0$ leads to
$$(2x^2-1)(16x^4-12x^2+1)=0\implies x=\pm \frac1{\sqrt2},\pm\frac{\sqrt5\pm 1}{4}$$
Thus, all the real solutions are:
$(0,0), (1,1), (-1,-1),
(\frac 1{\sqrt2},-\frac 1{\sqrt2}), (-\frac 1{\sqrt2},\frac 1{\sqrt2}), $
$(\frac {\sqrt5 + 1}{4},\frac {-\sqrt5 + 1}{4}), (\frac {\sqrt5 - 1}{4},\frac {-\sqrt5 - 1}{4}),
(\frac {-\sqrt5 + 1}{4},\frac {\sqrt5 + 1}{4}), (\frac {-\sqrt5 - 1}{4},\frac {\sqrt5 - 1}{4})$
Note that the system of the two equations are implicitly of ninth order. So, it is expected to have nine pairs of solutions.
This is NOT a full solution. Since contestants are asked to calculate the value of $p+q$ (without proof), they only need to obtain a numerical answer. I like to call this method "obtain answer by special case", where you can "tweak" the question to make it easier to calculate, as long as it still fits the constraints. If the question does indeed have a fixed solution, we will obtain the correct answer from our special case.
Since time is of the essence in competitions, this is used whenever possible (e.g. turning quadrilaterals into squares, letting constants be zero, etc.) I'd imagine the full solution to be quite a bit more complicated.
Best Answer
I guess I can explain their solution in more detail here*. So first (which they in fact do last in their solution), we want to get rid of the $2y+2x$ on the RHS of the second equation. We can do this via the slick substitution $(W,X,Y,Z)=(w-1,x+1,y+1,z-1)$ which gives us $$(W+1)(X+Y+Z-1) + (Y-1)(Z+X) + (X-1)(Z+1) = 2X-2+2Y-2+193$$ $$W(X+Y+Z)+XY+YZ+ZX-W-Z+X+Y-2=2X+2Y+193-4$$ And then we can add the first equation $W+X+Y+Z=25$ and we then get $$W(X+Y+Z)+XY+YZ+ZX=193+25-2=216$$ Now we are left with finding the maximum of $W$ with the above equation and the constraint $$X+Y+Z+W=25$$ We do this via the inequality $(X+Y+Z)^2\geq 3(XY+YZ+ZX)$, which is a direct result of the fact that $$(X-Y)^2+(Y-Z)^2+(Z-X)^2 \geq 0$$ This is all motivated because the term $X+Y+Z$ is quite nice since we can replace it by $25-W$. So now we have $$216 \leq W(25-W) + \frac{(25-W)^2}{3}$$ $$-75W+3W^2+216*3-25^2+50W-W^2 \leq0$$ $$2W^2-25W+23\leq0$$ Since the $W^2$ term has a positive coefficient, the $W$ which satisfy this inequality lie in the region between the two solutions, so the maximum solution is also the maximum $W$ and thus we have $$W\leq \frac{25+\sqrt{25^2-4*2*23}}{2*2}=\frac{25+\sqrt{441}}{4}=\frac{23}{2}$$ And then remembering to switch back to $w=W+1$ we have $$w\leq \frac{25}{2}$$ so our solution is $25+2=27$.
*(also, you changed their z to an x so I went along with that here)