$\newcommand{\angles}[1]{\langle #1 \rangle}$ Here is an idea which may motivate a real solution.

Claim. This is possible for $n = 2$ (i.e. in $\mathbb{R}^2$) with the added assumption that $\nabla f(0) + \nabla f(0)^T$ is a non-negative scalar matrix.

Define $T : \mathbb{R}^2 \to M_2(\mathbb{R})$ by \begin{equation}\tag{1} T(x) = \frac{1}{\lVert x \rVert^2} \begin{bmatrix} \angles{f(x), x} & -\angles{f(x), Rx} \\ \angles{f(x), Rx} & \angles{f(x), x} \end{bmatrix}, \end{equation} where $R : \mathbb{R}^2 \to \mathbb{R}^2$ is rotation by $\pi/2$. The motivation for this particular definition will be revealed later.

Of course, $T$ is not well-defined at $0$ – we really want $T$ to be the unique continuous extension to $\mathbb{R}^2$ of the function $\mathbb{R}^2 \setminus \{0\} \to M_2(\mathbb{R})$ defined by (1), and we must verify that this extension exists. Since the function is continuous on $\mathbb{R}^2 \setminus \{0\}$, it suffices to show that $$\lim_{x \to 0} \frac{1}{\lVert x \rVert^2} \begin{bmatrix} \angles{f(x), x} & -\angles{f(x), Rx} \\ \angles{f(x), Rx} & \angles{f(x), x} \end{bmatrix}$$ exists. To establish the existence of this limit, we must only verify that both $$\lim_{x \to 0} \frac{\angles{f(x), x}}{\lVert x \rVert^2}$$ and $$\lim_{x \to 0} \frac{\angles{f(x), Rx}}{\lVert x \rVert^2}$$ exist.

Now we note that $$\lim_{x \to 0} \left\lvert \frac{\angles{f(x) - \nabla f(0) x, x}}{\lVert x \rVert^2} \right\rvert \leq \lim_{x \to 0} \frac{\lVert f(x) - \nabla f(0) x \rVert \lVert x \rVert}{\lVert x \rVert^2} = \lim_{x \to 0} \frac{\lVert f(x) - \nabla f(0) x \rVert}{\lVert x \rVert} = 0$$ by Cauchy-Schwarz + the fact that $f$ is differentiable. Thus, $$\lim_{x \to 0} \frac{\angles{f(x) - \nabla f(0) x, x}}{\lVert x \rVert^2} = 0,$$ and for the same reason we have $$\lim_{x \to 0} \frac{\angles{f(x) - \nabla f(0) x, Rx}}{\lVert x \rVert^2} = 0.$$ Thus, $$\lim_{x \to 0} \frac{\angles{f(x), x}}{\lVert x \rVert^2} = \lim_{x \to 0} \frac{\angles{\nabla f(0) x, x}}{\lVert x \rVert^2} \qquad \text{and} \qquad \lim_{x \to 0} \frac{\angles{f(x), Rx}}{\lVert x \rVert^2} = \lim_{x \to 0} \frac{\angles{\nabla f(0) x, Rx}}{\lVert x \rVert^2}.$$

Here is where the assumption that $\nabla f(0) + \nabla f(0)^T$ is a scalar. Writing out what happens componentwise, we find that the first limit equals $\angles{\nabla f(0) e_1, e_1}$ and that the second limit equals $\angles{\nabla f(0) e_1, e_2}$, where $e_1 = (1,0)$ and $e_2 = (0,1)$.

Thus, $T : \mathbb{R}^2 \to M_2(\mathbb{R})$ is well-defined, with $$T(0) = \begin{bmatrix} \angles{\nabla f(0) e_1, e_1} & -\angles{\nabla f(0) e_1, e_2} \\ \angles{\nabla f(0) e_1, e_2} & \angles{\nabla f(0) e_1, e_1} \end{bmatrix}.$$ In fact, using the fact that $\nabla f(0) + \nabla f(0)^T$ is a scalar again, we have that $T(0) = \nabla f(0)$. In particular, $T(0)$ is positive semidefinite, because $T(0) + T(0)^T = \nabla f(0) + \nabla f(0)^T$ is PSD.

Moreover, $T(x)$ is PSD for all $x \neq 0$ because $$T(x) + T(x)^T = \frac{2}{\lVert x \rVert^2} \begin{bmatrix} \angles{f(x), x} & 0 \\ 0 & \angles{f(x), x} \end{bmatrix}$$ is a non-negative scalar matrix. Thus, $T$ takes values only in PSD matrices. Finally, we claim that $T(x)x = f(x)$. It suffices to verify this identity only when $x \neq 0$, since both sides are continuous in $x$. Here is where we (finally) motivate the definition of $T$:

Assuming $x \neq 0$, $\{x, Rx\}$ is a basis for $\mathbb{R}^2$. Thus, there is a unique matrix $M$ such that $Mx = f(x)$ and $M(Rx) = Rf(x)$. Writing down $M$ explicitly yields exactly the formula for $T(x)$ given in (1). In particular, we have $T(x)x = f(x)$ for all $x \neq 0$.

QED.

A few notes about this "non-negative scalar matrix condition":

1. It is always satisfied in $\mathbb{R}^1$ by the assumption that $f$ is differentiable and PSD.
2. It is satisfied in the example you gave with $f(x_1, x_2) = (-x_2^2, x_1 x_2)$ (since in this case $\nabla f(0) = 0$), and the formula (1) gives exactly the PSD matrix $\begin{bmatrix} 0 & x_2 \\ -x_2 & 0 \end{bmatrix}$ you suggested.

I'd like to conjecture that this is possible for $n = 2$, even without assuming the non-negative scalar matrix condition. Here is my reasoning:

1. For each $x \in \mathbb{R}^2$, there is at least one PSD matrix $T$ such that $Tx = f(x)$. Indeed, we may use (1) when $x \neq 0$ and choose $T = 0$ when $x = 0$.
2. For each $x \in \mathbb{R}^2$, let $\mathcal{Z}_x$ be the set of PSD matrices $Z \in M_2(\mathbb{R})$ such that $Zx = f(x)$. This set is always nonempty by the previous step.
3. $\mathcal{Z}_x$ "varies continuously in $x$", and thus there should be a continuous function $L : \mathbb{R}^2 \to M_2(\mathbb{R})$ such that $L(x) \in \mathcal{Z}_x$ for all $x$. To be more precise about the way that $\mathcal{Z}_x$ varies with $x$, note that $Zx = f(x)$ cuts out a continuously varying $2$-dimensional subspace of $M_2(\mathbb{R})$ (except when $x = 0$). To form the set $\mathcal{Z}_x$, we intersect this $2$-dimensional subspace with the set of PSD matrices. By Sylvester's Criterion, the set of PSD matrices is cut out by $3$ polynomial inequalities: a matrix $M$ is PSD iff all of the principal minors of $M + M^T$ are non-negative. That is, we must have $M_{1,1}, M_{2,2}, \det(M + M^T) \geq 0$. So overall, the set $\mathcal{Z}_x \subseteq M_2(\mathbb{R})$ is defined by finitely many polynomial inequalities of degree $\leq 2$, with the coefficients in these polynomials depending at most linearly on $x$. Since this behavior is relatively tame, it seems reasonable to me to claim that there must exist a continuous choice function, but I am not well-versed in this sort of thing. Maybe there's a pathology I'm overlooking!