This is a partial result towards an answer.
We can try a brute force approach to enumerate all functions which preserve conics by definition. We begin by recalling an ODE that all conics satisfy.
$$ C[y] = 9(y')^2 (y''''') -45(y'')(y''')(y'''') +40(y''')^3 = 0 $$
So given a parametric curve as $x(t)=t, y(t)$ as long as the $y(t)$ obeys the above equation then locally this curve is tracing out a conic.
Let us now consider some arbitrary function $\mathbb{R}^2 \rightarrow \mathbb{R}^2 : \phi = (\phi_x(x,y), \phi_y(x,y))$
This function has two components and accepts an $(x,y)$ pair and produces as $\phi_x, \phi_y$ output. If $\phi$ sends conics to conics then that means that if $(x(t), y(t))$ form a conic then $\phi_x(x(t), y(t)), \phi_y(x(t), y(t))$ must also form a conic. If we define a curve $(t, y(t))$ to form a conic IFF it obeys the equation above then the conclusion we have is that
$$ C[\phi_y] = C[y]*F(x,y,\phi_x ,\phi_y... ) $$
It's clear then that if $C[y] = 0$ then $C[\phi_y] = 0$, for any choice of $F$. So if $\phi_y$ obeys ANY such equation then it necessarily sends conics to conics.
Putting this together we have that
$$ 9(\phi_y(t, y(t))')^2 (\phi_y(t, y(t))''''') -45(\phi_y(t, y(t))'')(\phi_y(t, y(t))''')(\phi_y(t, y(t))'''') +40(\phi_y(t, y(t))''')^3 = \left( 9(y')^2 (y''''') -45(y'')(y''')(y'''') +40(y''')^3 \right) F(x,y,\phi_x, \phi_y, ... ) $$
Let's allow $F=1$ for now. We will revisit it later:
$$ 9(\phi_y(t, y(t))')^2 (\phi_y(t, y(t))''''') -45(\phi_y(t, y(t))'')(\phi_y(t, y(t))''')(\phi_y(t, y(t))'''') +40(\phi_y(t, y(t))''')^3 = \left( 9(y')^2 (y''''') -45(y'')(y''')(y'''') +40(y''')^3 \right) $$
Recall that these derivatives can be expanded via partial derivatives such as: $\phi_y(t, y(t))' = \partial_x[\phi_y] + \partial_y[\phi_y] y'(t) $
Thus the ENTIRE LHS can be expanded in terms of partial derivatives to look something as follows
$$ F_1 + F_2 y + F_3 y' + F_4 (y')^2 + ... $$
I.E. a polynomial in the derivatives of $y$ with coefficient $F_i$ being a function of partial derivatives of $\phi_y$.
This LHS has to equal $ \left( 9(y')^2 (y''''') -45(y'')(y''')(y'''') +40(y''')^3 \right) $ and the only way that can happen is if the $F_i$ corresponding to the $(y')^2 (y''''')$ term is equal to $9$, the $(y'')(y''')(y'''')$ term is equal to $-45$, the $(y''')^3$ is equal to $40$ and the rest are zero. So we end up with a massive system of partial differential equations almost ALL set to 0 except 3 lines:
$$ F_1 (\phi_y, \partial_x[\phi_y], \partial_y[\phi_y], ... ) = 0 \\
F_2 (\phi_y, \partial_x[\phi_y], \partial_y[\phi_y], ... ) = 0 \\
\vdots \\
F_{k_1} (\phi_y, \partial_x[\phi_y], \partial_y[\phi_y], ... ) = 9 \\
\vdots \\
F_{k_2} (\phi_y, \partial_x[\phi_y], \partial_y[\phi_y], ... ) = -45 \\
\vdots \\
F_{k_3} (\phi_y, \partial_x[\phi_y], \partial_y[\phi_y], ... ) = 45 \\
\vdots
$$
The solutions to THIS system are examples of functions which send conics to conics. Now some obervations:
As a system with WAY MORE equations than unknowns it is surprising ANY solutions exist at all.
This assumes $\phi_x(x,y) = x$ we know a larger solution space exists than this but we need a more generic $C$ operator to represent them
This is letting $F=1$. There are infinitely many other choices of $F$ that can be tried but they form a finite dimensional space in derivatives of $y$. If you pick a bad $F$ then this becomes unsolvable for obvious reasons.