$\color{green}{\textbf{UPDATED.}}$
Since $f'(x) \ge 0,$ then $f(x)$ is the non-decreasing function of $x.$
If $n=0,$ then the given inequality becames identity.
Since
$$\dfrac{2(xf(x)+yf(y))}{(x+y)(f(x)+f(y))}
= 1 + \dfrac{(x-y)(f(x)-f(y))}{(x+y)(f(x)+f(y))},\tag1$$
then the given inequality can be presented in the equivalent form of
$$f\left(\dfrac{x^2+y^2}{x+y}\right) \ge \left(1+\dfrac{(x-y)(f(x)-f(y))}{(x+y)(f(x)+f(y))}\right)f\left(\dfrac{x+y}2\right).\tag2$$
$\color{brown}{\textbf{Case 1.}\quad \mathbf{n>0,\quad a_0=0}.}$
$\textbf{Case 1.1.}\quad \mathbf{(x,y)\in [1,\infty)^2}.$
Denote
$$g(x) = \dfrac1x f(x)\ \Rightarrow f(x)=xg(x),$$
then
$$f(xy) = xy\, g(xy)\ge xy\,g(x) = y f(x),$$
$$2(x^2+y^2) = (x+y)^2+(x-y)^2,$$
$f\left(\dfrac{x^2+y^2}{x+y}\right) \ge \left(1+\dfrac{(x-y)^2}{(x+y)^2}\right)f\left(\dfrac{x+y}2\right).\tag3$
On the other hand, the inequality
$$1+\dfrac{(x-y)^2}{(x+y)^2}\ge1 + \dfrac{(x-y)(f(x)-f(y))}{(x+y)(f(x)+f(y))},\tag4$$
or
$$(x-y)^2(f(x)+f(y))\ge (x^2-y^2)(f(x)-f(y)),$$
$$2(x-y)(yf(x)-xf(y))\ge 0,$$
$$2xy(x-y)(g(x)-g(y))\ge 0,$$
is trivial for increasing $g(x).$
Multiplying of $(3)$ and $(4)$ leads to $(2).$
$\textbf{Case (1.1) is proved.}$
$\textbf{Case 1.2.}\quad \mathbf{m=\min(x,y) \in(0,1)}.$
Polynomial $f(x)$ can be presented in the form of
$$f(x) = g(mx),$$
where
$$g(u) = f\left(\dfrac um\right) = \sum\limits_{i=0}^n \dfrac{a_i}{m^i}\cdot u^i.$$
In accordance with the $\text{Case 1.1},\ \ g(u)$ satisfies to the given inequality if $(u,v)\in[1,\infty).$
$$(u+v)g\Big(\frac{u^2+v^2}{u+v}\Big)(g(u)+g(v))\geq 2(ug(u)+vg(v))g\Big(\frac{u+v}{2}\Big),$$
$$(u+v)f\Big(\frac 1m\frac{u^2+v^2}{u+v}\Big)
\left(f\Big(\frac um\Big)+f\Big(\frac vm\Big)\right)
\geq 2\left(uf\Big(\frac um\Big)+vf\Big(\frac um\Big)\right)f\Big(\frac{u+v}{2m}\Big),$$
$$\left(\Big(\frac um\Big)+\Big(\frac vm\Big)\right)
f\left(\frac{\Big(\frac um\Big)^2+\Big(\frac vm\Big)^2}{\Big(\frac um\Big)+\Big(\frac vm\Big)}\right)
\left(f\Big(\frac um\Big)+f\Big(\frac vm\Big)\right)\\
\geq 2\left(\Big(\frac um\Big)f\Big(\frac um\Big)
+v\Big(\frac vm\Big)\Big(\frac vm\Big)\right)
f\left(\frac{\Big(\frac um\Big)+\Big(\frac vm\Big)}{2}\right).$$
Therefore,
$$(x+y)f\Big(\frac{x^2+y^2}{x+y}\Big)(f(x)+f(y))\geq 2(xf(x)+yf(y))f\Big(\frac{x+y}{2}\Big),$$
where $(x,y) \in[m,\infty).$
Since $m$ is an arbitrary real number, then $\textbf{Cases (1.2) and (1) are proved.}$
$\color{brown}{\text{Case 2.}\quad \mathbf{n>0,\quad a_0>0}.}$
In this case polynomial $f(x)$ can be presented as the production of the factors $(x+r_j)$ and $((x+p_k)^2+q_k^2).$
Denote $c = \min(r_j,p_k).$
$\textbf{Case 2.1.}\quad\mathbf{c\in\{r_k\}}.$
The polynomial $f(x)$ can be presented in the form of
$$f(x) = (x+r)g(x+r),$$
wherein the polynomial $g(u)$ has positive coefficients.
If $r\in[1,\infty),$ then
$$f(xy) = (x+r)(y+r) g((x+r)(y+r))\ge (x+r)yg((x+r)(y+r) = yf(x)),\tag5$$
and $\textbf{the further proof is similar to the Case 1.1}.$
If $r\in[0,1)$ and $(x,y)\in[1-r,\infty)^2,$ then $(5)$ is correct,
and $\textbf{the further proof is similar to the Case 1.1}.$
If $r\in[0,1)$ and $(x,y)\in(0,1-r),$ then the polynomial $f(x)$ can be presented in the form of
$$f(x) = \dfrac{x+r}r g\left(\dfrac{x+r}r\right),$$
wherein the polynomial $g(u)$ has positive coefficients.
Then $\textbf{the further proof is similar to the Case 1.2}.$
$\textbf{Case 2.2.}\quad\mathbf{c\in\{p_j\}.}$
The polynomial $f(x)$ can be presented in the form of
$$f(x) = ((x+p)^2+q^2)g(x+p),$$
wherein the polynomial $g(u)$ have positive coefficients.
At the same time, $p$ can be negative.
For example, let us consider the polynomials in the form of
$$f^\,_k(x) = 1+x^{4k+2} = (1+x^2)(1-x^2+x^4-\dots+x^{4k}).$$
If $k=1,$ then $4k+2=6,$
$$x^4-x^2+1 = (x^2+1)^2-3x^2 = (x^2-\sqrt3\,x+1)(x^2+\sqrt3\,x+1)$$
$$=\left(x^2-2x\cos\dfrac\pi6+1\right)\left(x^2+2x\cos\dfrac\pi6+1\right)$$
$$=\left(\left(x-\cos\dfrac\pi6\right)^2+\sin^2\dfrac\pi6\right)
\left(\left(x+\cos\dfrac\pi6\right)^2+\sin^2\dfrac\pi6\right)$$
Then
$$f^\,_1(x)=(x^2-2zx+1)(1+x^2)(x^2+2zx+1) = (x^2-2zx+1)g^\,_1(x-z),$$
wherein $g^\,_1(u)$ is the polynomial with the positive coefficients.
Easy to show that for $k\ge1$
$$f^\,_k(x)=(x^2-2zx+1)g^\,_k(x-z),\tag6$$
where
$$z= \cos^2\dfrac\pi{4k+2}\tag7$$
and
$$g^\,_k(u)=\dfrac{1+(u+z)^{4k+2}}{1-z^2+u^2} = \prod\limits_{j=2}^k
\left(\left(u+z-\cos\dfrac{2j-1}{4k+2}\pi\right)^2+\sin^2 \dfrac{2j-1}{4k+2}\pi\right)$$
is the polynomial with the positive coefficients.
In particular, for
$$x=\dfrac z2(1+d),\quad y=\dfrac z2(1-d)$$
one can get
$$x+y = z,\quad x-y = zd,\quad x^2+y^2 = \dfrac{z^2}2(1+d^2),$$
$$f^\,_k\left(\dfrac z2(1+d^2)\right) - \left(1+d\dfrac{f^\,_k\left(\dfrac z2(1+d)\right)-f^\,_k\left(\dfrac z2(1-d)\right)}
{f^\,_k\left(\dfrac z2(1+d)\right)+f^\,_k\left(\dfrac z2(1-d)\right)}\right)f^\,_k\left(\dfrac z2\right)$$
$$=\left(1+\left(\dfrac z2\right)^{4k+2}(1+d^2)^{4k+2}\right)
- \left(1+d\cdot\dfrac{\left(\dfrac z2\right)^{4k+2}\left((1+d)^{4k+2}-(1-d)^{4k+2}\right)}
{2+\left(\dfrac z2\right)^{4k+2}\left({\large\mathstrut}(1+d)^{4k+2}+(1-d)^{4k+2}\right)}\right)$$
$$=\left(\dfrac z2\right)^{4k+2}(1+d^2)^{4k+2}\left(1-\dfrac{d\left((1+d)^{4k+2}-(1-d)^{4k+2}\right)}
{(1+d^2)^{4k+2}\left(2+\left(\dfrac z2\right)^{4k+2}\left({\large\mathstrut}(1+d)^{4k+2}+(1-d)^{4k+2}\right)\right)}\right)$$
$$\le \left(\dfrac z2\right)^{4k+2}(1+d^2)^{4k+2}\left(1-\dfrac{d\,(1+d)^{4k+2}}
{(1+d^2)^{4k+2}\left(2+\left(\dfrac{1+d}2z\right)^{4k+2}\right)}\right),$$
and this allows to find counterexamples
- $k=2,\quad 4k+2=10,\quad (x,y) = \dfrac z2(1\pm d),\quad$
where$\quad z = \cos \dfrac\pi{10}\approx0.951,\quad d\in[0.328, 0.774]\quad$
(see also Wolfram Alpha calculations)

- $k=3,\quad 4k+2=14,\quad (x,y) = \dfrac z2(1\pm d),\quad$
where$\quad z = \cos \dfrac\pi{14}\approx0.975,\quad d\in[0.229, 0.856]\quad$
(see also Wolfram Alpha calculations)

- $k=4,\quad 4k+2=18,\quad (x,y) = \dfrac z2(1\pm d),\quad$
where$\quad z = \cos \dfrac\pi{10}\approx0.985,\quad d\in[0.181, 0.892]\quad$
(see also Wolfram Alpha calculations)

Therefore, in the case $(2.2)$ the given inequality $\color{brown}{\textbf{ allows counterexamples}}.$