An inequality for polynomials with positives coefficients
Even for $f(x)=x^n$ it's not so easy.
We need to prove $$(x+y)\Big(\frac{x^2+y^2}{x+y}\Big)^n(x^n+y^n)\geq 2\left(x^{n+1}+y^{n+1}\right)\Big(\frac{x+y}{2}\Big)^n$$ or $$2^{n-1}(x^2+y^2)^n(x^n+y^n)\geq\left(x^{n+1}+y^{n+1}\right)(x+y)^{2n-1}.$$ Now, let $x=ty$.
Also, since our inequality is symmetric, we can assume that $t\geq1.$
Thus, we need to prove that $g(t)\geq0,$ where $$g(t)=(n-1)\ln2+n\ln(t^2+1)+\ln(t^n+1)-\ln\left(t^{n+1}+1\right)-(2n-1)\ln(t+1).$$ Now, $$g'(t)=\frac{h(t)}{(t^2+1)\left(t^{n+1}+1\right)\left(t^n+1\right)(t+1)},$$ where $$h(t)=n(t-1)^3(t+1)t^{n-1}+2n(t-1)\left(t^{2n+1}+1\right)-(t^2+1)\left(t^{2n}-1\right).$$ Now, prove that $$h(1)=h'(1)=h''(1)=0$$ and $h'''(t)\geq0$ for all $t\geq1.$
Here is a partial answer. I show below that your inequality holds when $n\leq 5$.
Let $d=(x+y)^n\frac{LHS-RHS}{(y-x)^4}$. This turns out to be a polynomial, and a computation in PARI-GP (see below) reveals that $d$ is a polynomial in positive coefficients in $x,y,a_0,a_1,\ldots,a_n$ when $n\leq 4$ ; I call $d$ a completely positive polynomial.
Things become more complicated for $n=5$, because we have a "negative component" equal to
$$ -\frac{35}{16}a_0a_5(x^4y^3+x^3y^4) $$
However, we are still able to obtain positivity of $d$ in this case because the $a_0a_5$ coefficient is of the form $(y-x)^2$ times something completely positive.
As $n$ grows, the negative monomials grow in number too but they remain a minority, so the method for $n=5$ can probably be generalized.
Here is the PARI-GP code I used :
n0=4
aa(k)=eval(Str("a",k))
expr1=sum(j=0,n0,aa(j)*(u^j))
expr2=subst(expr1,u,p/q)*(q^n0)
expr3=subst(subst(expr2,p,x^2+y^2),q,x+y)
exprx=subst(expr1,u,x)
expry=subst(expr1,u,y)
main0=(x+y)*expr3*(exprx+expry)-\
2*((x+y)^n0)*(x*exprx+y*expry)*(subst(expr1,u,(x+y)/2))
main=main0/((y-x)^4)
I have tried many approaches to prove this inequality, but none worked. By now, I think that the inequality doesn't hold and thus I started to search for a counterexample. Michael Rozenberg has gave a proof for the special case $f(x) = x^n$ and Ewan Delanoy verifies this inequality for polynoms with degree at most $5$.
First we note that the condition that all $a_i >0$ are positive is unnecessary: If the inequality is valid for all $a_i> 0$, then we could let $a_i \downarrow 0$ for any index $i =0,\ldots,n$ and the inequality would remain to be valid.
However, it is in general wrong: Let $x=1$ and $y=t$ and take $f(x) = 1+x^{10}$. Then the function $$g(t):= (1+t) f \Big( \frac{t^2+1}{1+t} \Big)(f(1)+f(t)) - 2 (f(1)+t f(t)) f \Big( \frac{t+1}{2} \Big) $$ is negative for $t=0.5$: See this plot in WolframAlpha.
Conjecture: Is the inequality valid if we additionally require that $a_0 \ge a_1 \ge ... \ge a_n$?
I couldn't find any counterexample on this strengthened variant, but also no promising approach to prove this. Maybe someone has an idea?