Find all polynomials $P, Q$ such that $(x^2+ax+b)P(x^2+cx+d)=(x^2+cx+d)Q(x^2+ax+b)$ 
Find all polynomials $P(x), Q(x)$ such that
$(x^2+ax+b)P(x^2+cx+d)=(x^2+cx+d)Q(x^2+ax+b)$
where $a, b, c, d$ are all different real numbers.

I found this problem while solving functional equations in polynomials. Here's what I've found so far:
Suppose $u$ is a root of $x^2+ax+b$ which is not a root of $x^2+cx+d$. Then, substituting $x=u$ gives $Q(0)=0$.
Similarly, suppose $v$ is a root of $x^2+cx+d$ which is not a root of $x^2+ax+b$. Then, substituting $x=v$ gives $P(0)=0$.
Therefore, $0$ is a root of both $P$ and $Q$, so we can write $P(x)=xf(x)$ and $Q(x)=xg(x)$ for some polynomials $f, g$.
The original equation then becomes
$h(x)f(x^2+cx+d)=h(x)g(x^2+ax+b)$, where
$h(x)=(x^2+ax+b)(x^2+cx+d)$.
After canceling out $h(x)$, we get
$f(x^2+cx+d)=g(x^2+ax+b)$.
However, I can't find a way to make further progress. What should I do to finish solving this problem?

As Ng Chung Tak said in the comments above, the only solutions are $P = Q = \lambda X$, $\lambda\in \mathbb{R}$, if $a \neq c$ (though the answer seems more complicated for $a = c$, but this is not the object of this question).
Let $(P,Q) \in \mathbb{R}[X]^2$ be a solution of our polynomial equation, equation which we'll name $(*)$.
By taking the degrees in $(*)$, we get:
$$2 + 2\deg(P) = 2 + 2\deg(Q) \quad\Rightarrow\quad \deg(P) = \deg(Q) =: n$$
Similarly, by looking at the dominant coefficients of the polynomials on each side of $(*)$, we also have that $P$ and $Q$ have the same leading coefficient $\lambda$.
If $n = -\infty$, and thus $P = Q = 0$, we are done, and it is easily shown that $n$ can't be $0$ otherwise $a = c$ and $b = d$, so we can suppose $n \geq 1$.
Let $P =: \sum_{k = 0}^n \alpha_k X^k$ and $Q =: \sum_{k = 0}^n \beta_k X^k$ with $\alpha_n = \beta_n = \lambda \neq 0$. We then obtain:
$$\begin{align}(X^2 + aX + b)P(X^2 + cX + d) &= (X^2 + aX + b)\sum_{k=0}^n \alpha_k (X^2 + cX + d)^k =: S \\ (X^2 + cX + d)Q(X^2 + aX + b) &= (X^2 + cX + d)\sum_{k=0}^n \beta_k (X^2 + aX + b)^k =: T\end{align}$$
Let's look at the coefficient of $S$ that's in front of $X^{2n+1}$.
The highest power of $X$ appearing in $(X^2 + cX + d)^k$ for $k \leq n - 1$ is $X^{2n-2}$, which multiplied by $X^2$ only grants $X^{2n}$, hence only $(X^2 + cX + d)^n$ contributes, and the coefficient of $S$ in front of $X^{2n+1}$ is thus $\lambda (nc + a)$ (the only way to get $X^{2n+1}$ is to use $n$ times a factor in $X^2$ and exactly once a factor in $X$). In the same fashion, the corresponding coefficient for $T$ is $\lambda(na + c)$.
Therefore, since $(*)$ is just $S = T$:
$$\lambda(nc+a) = \lambda(na+c) \quad\Rightarrow\quad nc+a = na+c \quad\Rightarrow\quad n(c-a) = c - a \quad\Rightarrow\quad n = 1$$
As such, we have: $P = \lambda X + \alpha_0$ and $Q = \lambda X + \beta_0$, and so:
$$\begin{align}S &= (X^2 + aX + b) \big(\lambda(X^2 + cX + d) + \alpha_0\big) \\ T &= (X^2 + cX + d) \big(\lambda(X^2 + aX + b) + \beta_0\big)\end{align}$$
Which means that, by simplifying by $\lambda (X^2 + aX + b)(X^2 + cX + d)$:
$$\alpha_0(X^2 + aX + b) = \beta_0(X^2 + cX + d)$$
We can conclude that $\alpha_0 = \beta_0$ by looking at the coefficients in front of $X^2$, and then that $\alpha_0 a = \alpha_0 c$ by looking at $X$, which finally provides: $\alpha_0 = 0$, and so: $P = Q = \lambda X$.
Reciprocally, those are indeed solutions, and they are therefore exactly the solutions to $(*)$ when $a \neq c$ (not even any need to know anything about $b$ and $d$).

The case $a = c, b \neq d$ (which I'm aware is not really part of this question, but I still want to talk briefly about it) is actually way less trivial it seems:
Take $a = c = d = 0, b = -1$. You can prove that for any $\alpha \in \mathbb{R}$, the couple $(X^2 + \alpha X, X^2 + (\alpha + 1)X)$ is solution of $(*)$ (there's not even uniqueness modulo a multiplicative constant per degree in this case, $P$ is different from $Q$, and so on...). I wonder if solutions like this exist for all quadruplets such that $a = c$?
Moreover, the proof for $a \neq c$ can be generalised to all integral domains of characteristic zero, which $\mathbb{R}$ and $\mathbb{C}$ are, but not so much to rings of positive characteristic, though the common degree $n$ can at least be reduced to $1$ modulo said characteristic.
