It's not hard to see that $H$ must be of the form $\alpha ux + \beta uz + \gamma yx + \delta yz$ by $\mathbb{R}$-linearity of $F$ (see Jan-Cristoph Schlage-Puchta's answer for a fuller explanation).
By using symmetry of multiplication, we can assume without loss of generality that $\beta = \gamma$.
We have the equation $$F(P \times Q) = H(F(P), F(Q), P, Q) = \alpha F(P) F(Q) + \beta F(P) Q + \gamma P F(Q) + \delta PQ$$
Setting $P = Q = 1$, we get $F(1) = \alpha F(1)^2 + \beta F(1) + \gamma F(1) + \delta$. Degree-checking tells us that if $F(1)$ isn't constant, then $\alpha = 0$. More careful checking tells us that $\beta + \gamma = 1, \delta = 0$. This leads to the solution $F_{M, R}(P) = P \times R$, where $R = F(1)$. That answer aside, we assume $F(1)$ is constant.
I'm going to take a detour for a moment, to talk about an action of the two-dimensional nonabelian Lie group on the space of pairs $(F, H)$ that satisfy $F(P \times Q) = H(F(P), F(Q), P, Q)$. This action comes from the map $g_{(a, b)} F = aF + b \operatorname{Id}$. We want the equation to remain the same, which leads us to:
$$(gH)(u, x, y, z) = aH(\frac{u - by}{a}, \frac{x - bz}{a}, y, z) + b yz$$
$$ = \frac{\alpha}{a} ux + (\beta - \frac{b}{a} \alpha) uz + (\gamma - \frac{b}{a} \alpha) yx + (a \delta - b \beta - b \gamma + \frac{b^2}{a} \alpha + b) yz$$
Back from the detour: we now set just $Q = 1$, getting $F(P) = \alpha F(P) F(1) + \beta F(P) + \gamma P F(1) + \delta P$. Rearranging, we get $(1 - \alpha F(1) - \beta) F(P) = (\gamma F(1) + \delta) P$. The obvious solutions have $F(P) = \lambda P$ for some $\lambda \in \mathbb{R}$, which have been discussed in other answers; otherwise, we must have $1 - \alpha F(1) - \beta = \gamma F(1) + \delta = 0$.
We now split into two cases: either $\alpha = 0$, or $\alpha \neq 0$. If $\alpha = 0$, then $\beta = \gamma = 1$. Using the action of the group, we can reduce to the case that $\delta = 0$, implying that $F(P \times Q) = F(P) Q + P F(Q)$, i.e. that $F$ is a derivation. This can only happen when $F(P) = R \times \partial P$ for some polynomial $R$. Undoing the group action, we get $F_{D, R, \lambda}(P) = \lambda P + R \times \partial P$.
If $\alpha \neq 0$, then we can use the action of the group to reduce to the case $\alpha = 1, \beta = \gamma = 0$. But then $\delta = 0$, giving us the equation $F(P \times Q) = F(P) F(Q)$ - in other words, $F$ is a homomorphism. Homomorphisms from $\mathbb{R}[X]$ are compositions: $F(P) = P \circ R$ for some polynomial $R$. Undoing the group action, we get $F_{H, R, \lambda} = \lambda P + P \circ R$.
So all solutions $(F, H)$ are of the forms:
$F(P) = \lambda P, \frac{1}{\lambda} \alpha + \beta + \gamma + \lambda \delta = 1$
$F(P) = \lambda P + R \times \partial P, \alpha = 0, \beta = \gamma = 1, \delta = \lambda$
$F(P) = \lambda P + c (P \circ R)$ with coefficients that aren't difficult to determine using the group action.