Here are some solutions.
Least Algebraic:
Let the foci of the parabolas be F and G. Let A,B be two of the points of intersection and C the midpoint of AB. Drop perpendiculars from C to the directrices of the two parabolas intersecting the curves at P and Q and the directrices at R and S. Then the tangents to the two parabolas at P and Q are parallel to AB. Draw perpendiculars to the tangents at P and Q that meet the axes of their respective parabolas at L and T, and draw the perpendicular to AB at C so that it meets the two axes in M and N as in the figure.
Then GLQS, GMCQ, FRPT, FRCN are parallelograms. Let O be on MN such that the distance to the axis FN is HO=UG. Then NO=CM so PMOT is also a parallelogram, so the distance from O to the axis GM is OK=FV. The point O so defined is on the perpendicular bisector of the segment joining each pair of intersection points, so they are chords of a common circle.
Medium algebraic:
This is essentially @Day Late Don's solution with a picture. Let the distance from the directrix of one parabola to its focus F be $h$ and the distance to the axis of the other parabola be $g$, and the respective distances for the other parabola with focus G be $j$ and $f$. Let A be one point of intersection and the distances from A to the axes be $x$ and $y$ as in the diagram. Let O be the point distance $f+j$ from the the directrix of the parabola with focus G and distance $h+g$ from the other directrix.
Then applying the Pythagorean theorem to right triangles GPA, FQA and AKO we have
$$
\begin{align}
(x+f-j)^2+y^2 & =(x+f)^2\\
y^2+j^2 & = 2j(x+f) \\
(g-y-h)^2+x^2 & = (g-y)^2\\
x^2+h^2 & = 2h(g-y) \\
AO^2 & = (x-j)^2+(y+h)^2\\
& = x^2+y^2+j^2+h^2-2xj+2yh \\
& = 2(jf+gh)
\end{align}
$$
The last expression does not depend on A, so the distance from each point of intersection to O is the same, thus they are concyclic. This looks like a lot of algebra, but I included it because it only really needs the Pythagorean theorem and sums and differences, so in theory could be eked out with a straight edge and compass.
More algebraic:
I thought this looked promising but is missing a non-algebraic proof that centroid of the quadrilateral defined by the intersection points is at the intersection of the axes. This follows easily from the algebraic approach as in the question (axes are $y=0$ and $x=0$), since the quartics have zero cubic terms, so the sum of the roots is $0$. I haven't found a more-graphical proof of this step.
Nevertheless, given that, let AB and CD be segments joining pairs of intersection points and M and N their respective midpoints, so that midpoint of MN is the centroid at the intersection of the axes. Drop perpendiculars from M and N to the directrix of one of the parabolas meeting the curve at Q and P.
Then AB is parallel to the tangent at Q and CD is parallel to the tangent at P. Since M and N are equidistant from each axis, P and Q are reflections in the axis of the parabola they're on, as are the tangents at P and Q.
Thus $\angle ABX = \angle DCW$ and a similar argument the other way gives $\angle YCB = \angle DAX$, and it follows that $\angle DAB + \angle BCD = \pi$ and ABCD is cyclic.
The following example is a supplement to achille hui's definitive answer.
Consider the family of ellipsoids
$${\cal E}_\lambda:\quad x^2-2\lambda xy+y^2+(1-\lambda^2)z^2=1-\lambda^2\qquad(-1<\lambda<1)\ .\tag{1}$$
Note that ${\cal E}_0$ is just the unit sphere.
We now project ${\cal E}_\lambda$ along the $y$-axis onto the $(x,z)$-plane. The points on ${\cal E}_\lambda$ giving rise to a boundary point of ${\cal Y}_\lambda:=\pi_y({\cal E}_\lambda)$ are the points where
$${\partial\over\partial y}\bigl(x^2-2\lambda xy+y^2+(1-\lambda^2)z^2\bigr)=2y-2\lambda x=0\ .\tag{2}$$
The points on ${\cal E}_\lambda$ satisfying $(2)$ are the points $(x,\lambda x,z)$ that also satisfy $(1)$, which amounts to the condition $$x^2+z^2=1\ .$$
It follows that ${\cal Y}_\lambda$ is the unit disk in the $(x,z)$-plane, irrespective of the value of $\lambda$.
By symmetry it follows that the projected ellipsoid ${\cal X}_\lambda:=\pi_x({\cal E}_\lambda)$ is the unit disk in the $(y,z)$-plane, irrespective of the value of $\lambda$.
Best Answer
I'd like to fill in Jack's remark that this is simply an eigenvector problem, though I won't do so from a variational perspective. Let's write the ellipse in bilinear form using matrices:
\begin{align} \mathbf{x}^T M \mathbf{x} &=\begin{pmatrix} x&y&z \end{pmatrix}\begin{pmatrix}\alpha_1&\beta_3&\beta_3\\ \beta_{2}&\alpha_2 & \beta_1\\ \beta_2 &\beta_1 & \alpha_3\end{pmatrix}\begin{pmatrix}{x\\y\\z}\end{pmatrix}\\\\ &=\alpha_1 x^2+\alpha_2 y^2+\alpha_3 z^2+2\beta_1 zy+2\beta_2 xz+2\beta_3 xy=1 \end{align}
Suppose the matrix $M$ has full rank i.e. three orthonormal eigenvectors $\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3$ with eigenvalues $\lambda_1,\lambda_2,\lambda_3$ respectively. Then $M$ may be diagonalized using the matrix $S=(\mathbf{v}_1\,\mathbf{v}_2\,\mathbf{v}_3)$ as
\begin{align} S^TMS &=\begin{pmatrix}\mathbf{v}_1^T\\\mathbf{v}_2^T\\\mathbf{v}_3^T\end{pmatrix} \begin{pmatrix}M \mathbf{v}_1&M\mathbf{v}_2&M \mathbf{v}_3\end{pmatrix}\\ &=\begin{pmatrix}\mathbf{v}_1^T\\\mathbf{v}_2^T\\\mathbf{v}_3^T\end{pmatrix} \begin{pmatrix}\lambda_1 \mathbf{v}_1&\lambda_2 \mathbf{v}_2 &\lambda_3 \mathbf{v}_3 \end{pmatrix}=\begin{pmatrix} \lambda_1&0&0\\0&\lambda_2&0\\0&0&\lambda_3\end{pmatrix} \end{align} where the last equality follows from orthonormality. Hence $S^T MS=D$ where $D$ is diagonal. But $S^T S=I$ (since $(S^T S)_{ij}=\mathbf{v}_i^T \mathbf{v}_j=\delta_{ij})$ by orthonormality), so we may write $M=S D S^T$ and therefore the ellipse becomes $$\mathbf{x}^T M \mathbf{x}=\mathbf{x}^T (S D S^T)\mathbf{x}=\mathbf{X}^T D \mathbf{X}=\lambda_1 X^2+\lambda_2 Y^2+\lambda_3 Z^2=1$$ where $\mathbf{X}=S^T \mathbf{x}=(\mathbf{v}^T_1 \mathbf{x},\mathbf{v}^T_2 \mathbf{x},\mathbf{v}^T_2 \mathbf{x})^T=(X,Y,Z)^T$.
In these coordinates, we just have an ellipsoid with principal $XYZ$-axes; in other words, its axes are lines whose tangent vectors are just the basis vectors $\hat{e}_X,\hat{e}_Y,\hat{e}_Z$. We can map this back to the original coordinates via $\mathbf{x}=S\mathbf{X}$. But $S\cdot (\hat{e}_X,\hat{e}_Y,\hat{e}_Z)=SI=S$, so the column vectors of $S$ are the images of the basis $XYZ$-vectors.
So we arrive at a cute conclusion: The tangent vectors of the principal axes of the ellipse are the eigenvectors of the coordinate matrix of the ellipse. Glancing at the $XYZ$-form of the ellipse, we conclude that we can also read off the length of these axes from the cooresponding eigenvalues. So the analytic-geometry problem of describing the ellipse is indeed equivalent to an eigenvalue-eigenvector problem in linear algebra.