Continuing the epic quest on finding matrix representations from here: Representation of hyperbolas. with a last part, the only conic section left: the parabola.
I will present one idea of how to structure this representation. But I hope to get examples of more elegant representations, hopefully ones which are conceptually similar to those of ellipses and hyperbolas.
Just to clarify: the objective is to find a matrix representation for a generating element such that the coordinates of a parabola are directly readable from some of the elements of a matrix power.
Remember addition can be represented as the block ${\bf M_{a+}} = \left(\begin{array}{cc} 1 & 0\\ a & 1 \end{array}\right)$ as we see that: $${\bf M_{a+}M_{b+}} = \left(\begin{array}{cc} 1 & 0\\ a & 1 \end{array}\right)\left(\begin{array}{cc} 1 & 0\\ b & 1 \end{array}\right) = \left(\begin{array}{cc} 1 & 0\\ a+b & 1 \end{array}\right)$$
Then observe that $({\bf M_{\Delta_t+}})^n$ will store $n\Delta_t$ at bottom left. And then expand the building block above with the row: $(\Delta_t\,1\,1)$, and we choose the column so we preserve triangularity.
If we take a look at this new matrix: $$\left(\begin{array}{ccc}1&0&0\\\Delta_t&1&0\\\Delta_t&1&1\end{array}\right)$$ Now the lower right 2x2 block will produce $n$ at position (3,2) if the matrix is raised to $n$.
One should be able to show that (3,1) becomes $\Delta_t \sum_{i=1}^{k} i$ which is the famous arithmetic sum which is a quadratic function. Now we have access to constant values, linear values and squares in our matrix.
What remains is to linearly combine these. We therefore create one last row and there we place the coefficients for our second degree polynomial $c+bx+ax^2$ : $(c+\Delta_t^2, b+\Delta_t, a+2\Delta_t)$. Our polynomial value will now end up in (4,1) such that $({\bf M}^n)_{(4,1)} = a(n\Delta_t)^2+b(n\Delta_t)+c$
All in all the generating element matrix looks like this: $${\bf M} = \left(\begin{array}{cccc}1&0&0&0\\\Delta_t&1&0&0\\\Delta_t&1&1&0\\\Delta_t^2+c&(\Delta_t+1)^2b&2a\Delta_t&0\end{array}\right)$$
We calculate the matrix powers, extract and plot $$x(n\Delta_t) = ({\bf M}^n)_{(2,1)}\\ y(n\Delta_t) = ({\bf M}^n)_{(4,1)}$$ which are compared to a "ground-truth" polynomial.
As you can see this is a quite big and cumbersome matrix as compared to the 2x2 generating elements for ellipses and hyperbola in the previous question. Does anyone have any idea of a more compact matrix representation?

As can be seen in second image there is also some minor flaw, likely a small constant fault.
If you make your conditions a bit different, it gets mathematically a bit more uniform. Instead of requiring $M^n$ having x,y components in it, you can instead require that if $\vec{r}$ is on the curve, then $M\vec{r}$ is also on the curve (and therefore $M^n\vec{r}$ too). A further simplification is to consider infinitesimal moves.
For a circle, that's easy. You start with $x^2+y^2=C$ (remember, if you start on a non-unit circle, you will stay on a non-unit circle, you describe the whole family with this). Then differentiating will get you to
$$x\,dx+y\,dy=0$$ Because you need a linear map without a constant factor (multiplication by $M$), you require $$(dx,dy)=\sigma\cdot (x,y)dt$$ leading to $$(x,y)\sigma(x,y)=0$$ A skew symmetric $M$ solves this nicely, $\sigma=\begin{bmatrix}0&-1\\1&0\end{bmatrix}$: $$(x,y)\to (x,y)+\sigma(x,y)dt=\underbrace{(1+\sigma dt)}_M(x,y)$$ In the limit $dt\to 0$, $M^n$ are rotational matrices with sines and cosines ($\sigma$ is the generator of Lie algebra). For symmetric hyperbolae, exactly the same procedure returns $\sigma=\begin{bmatrix}0&1 \\ 1&0\end{bmatrix}$. For ellipses and general hyperbolae, you the $1$ have to be replaced by the squared semiaxes.
Parabolae cannot be described as quadratic forms in 2D (they are nonhomogeneous), but you can go into the projective space which requires 3 components. In that space, you deal with $(x,y,z\equiv 1)$, and the condition is
$$ax^2+bxz+cz^2-zy=0$$ Differentiation: $$2ax\,dx+bx\,dz+bz\,dx+2cz\,dz - z\,dy-y\,dz=0$$ $$(x,y,z)\underbrace{\begin{bmatrix}2a & 0 & b \\ 0& 0& -1 \\ b & -1 & 2c\end{bmatrix}}_A(dx,dy,dz)=0$$ Again, we need a matrix $\sigma$ with the property $$(dx,dy,dz)=\sigma (x,y,z)$$ Now the condition is again that $B=A\sigma$ is antisymmetric. The choice of it is however not arbitrary, because we require $dz=0$ to stay on the same projection plane. This leads to
$$\sigma=\begin{bmatrix} 0 & 0 & 1 \\ 2a & 0 & b \\ 0 &0 & 0\end{bmatrix}$$
$$ M = 1+\sigma dt$$
You can't really say it's surprising, because that just means $$x\to x+dt$$ $$y\to y+(2ax+b)dt$$ We could have written this down immediately from the derivative (it's the obvious parameterization $x=t$, but I wanted to tell the procedure first.
Note that this is exact only in $dt\to 0$ limit. For a finite $\Delta t$, we get $$M^n=\begin{bmatrix}1&0&n\Delta t\\2an \Delta t&1&b n \Delta t +an(n+1) \Delta t^2\\0&0&1\end{bmatrix}$$ or written out, $$x = x_0+n \Delta t$$ $$y = y_0+(2a x_0+b) n \Delta t + a n(n+1) \Delta t^2$$ which satisfies a slightly modified quadratic equation. You can go back and change the $b \to b-a \Delta t$ coefficient of $\sigma$ to make it right in the finite $\Delta t$ case.
In the limit $n\to \infty$, $n\Delta t=t$, we get the continuous representation $$\begin{bmatrix}1&0&t\\2at & 1 & at^2+bt \\ 0 & 0 & 0\end{bmatrix}$$
To sum up: start with a $(x,y,1)$ vector on your curve. Multiply by $M$ to move to different points on the parabola. Every conic section can be handled this way.
Double-check my procedure, it is possible that I made some mistakes.