Quadratic P.S.D. differential operator that is invariant under $\textrm{SL}(2, \mathbb{R})$
Solution 1:
I will take your question seriously but not literally since you are asking for a differential operator $D: L^2({\mathbb R}^2)\to L^2({\mathbb R}^2)$ and there are no differential operators $D$ of order $>0$ which take all $L^2$-functions to $L^2$-functions: You would need distributions as values of $D$.
Thus, I will assume that $L^2({\mathbb R}^2)$ in your question means $C^\infty({\mathbb R}^2)$. Then you get your example:
Take $Z=x\frac{\partial}{\partial x} + y \frac{\partial}{\partial y}$: As a vector-field, it sends each point with coordinates $(x,y)$ to vector with coordinates $(x,y)$, this is why $Z$ is invariant under the action of $GL(2, {\mathbb R})$. Then take $D=Z\otimes Z$. As a differential operator, it acts on smooth functions by
$$
D: f\mapsto (x\frac{\partial f}{\partial x} + y \frac{\partial f}{\partial y})^2.
$$
One can prove that among strictly 1st order strictly quadratic PSD differential operators, up to scalar, this is the only one which is $SL(2, {\mathbb R})$-invariant.
The same works in higher dimensions as well, your $GL(n, {\mathbb R})$-invariant differential operator will be $$ Z\otimes Z, Z=\sum_{i=1}^n x_i\frac{\partial}{\partial x_i}. $$
If I were to take your question literally but not seriously, my answer would be $$ D: L^2({\mathbb R}^2)\to L^2({\mathbb R}^2), D(f)=a f^2 $$ where $a\ge 0$ is a fixed constant. Such $D$ is a PSD, continuous, quadratic differential operator of order 0.