The condition on your coefficients is an underdetermined system of quadratic equations. So specializing some of them (for instance setting them $0$), you can try to get a system with a finite number of (complex) solutions, and the chances are good that they are actually rational. Such systems can be handled with SageMath (which uses Singular, msolve or other systems as backend). In the specific case, a possible solution is \begin{align} Q_0 &= X^2 - XY + Y^2 + XZ\\\\ Q_1 &= X^2 - XY + YZ\\\\ Q_2 &= X^2 - XZ + YZ + Z^2. \end{align} Added later upon request in the comment: Here is the rather naive and straightforward Sage code which runs 160 seconds on my machine. It shows that a rational solution requires at least $11$ nonzero coefficients.
e2 = [(i, j) for i in range(3) for j in range(3-i)]
n2 = len(e2)
R = PolynomialRing(QQ, 'a', 3*n2)
Rxy.<x, y> = R[]
T0, T1, T2 = [sum(R.gen(m+k*n2)*x^i*y^j for m, (i, j) in enumerate(e2))
for k in range(3)]
f = x^3*y + y^3 + x
l = (T0*T2-T1^2-f).coefficients()
N = R.ngens()
for m in range(1, N):
print('m = ', m)
for s in Subsets(range(N), m):
i = min(s)
if i >= n2:
continue
l0 = l + [R.gen(i)-1] + [R.gen(j) for j in range(N) if
j != i and not j in s]
I = ideal(l0)
if I.dimension() == 0:
V = I.variety()
if len(V) > 0:
break
else:
continue
break
v = V[0]
Q0, Q1, Q2 = [sum(c.subs(v)*xy for c, xy in Q) for Q in [T0, T1, T2]]
print(f'Q0 = {Q0}\nQ1 = {Q1}\nQ2 = {Q2}\n{Q0*Q2-Q1^2 == f}')