I conjecture a stronger statement: that the function $F(x,y)$ of two variables $x,y$ obtained from $F_n(x)$ by substituting $y$ for $x^{2n}$ is convex for $0<y<x<0$. While this reduces to checking positive definiteness of a matrix of bivariate polynomials of degree about 20, obtained from the Hessian of $F$ by clearing common (positive) denominator, this should be possible to make to work by standard real algebraic geometry tools. Here is an (ugly) plot of $F$:

obtained by exporting the plot from jmol in Sage(math)

sage: var('x y')
sage: F(x,y)=log((x+y^2)*(1+y)*(1-x*y)/(x*(1+x*y)*(1-x^2*y)))
sage: plot3d(F,[0,1],[0,1])

EDIT: comments say that as $y\to 0+$, the Hessian becomes indefinite. That is, for my original claim to hold, one need to assume $y>\epsilon>0$, with $\epsilon$ possibly depending upon $x$.

I conjecture a stronger statement: that the function $F(x,y)$ of two variables $x,y$ obtained from $F_n(x)$ by substituting $y$ for $x^{2n}$ is convex for $0<y<x<0$.

EDIT2: this conjecture is wrong as stated. See e.g. the plot in the answer by Yaakov Baruch, or rotate the plot in Sage...

While this reduces to checking positive definiteness of a matrix of bivariate polynomials of degree about 20, obtained from the Hessian of $F$ by clearing common (positive) denominator, this should be possible to make to work by standard real algebraic geometry tools. Here is an (ugly) plot of $F$:

obtained by exporting the plot from jmol in Sage(math)

sage: var('x y')
sage: F(x,y)=log((x+y^2)*(1+y)*(1-x*y)/(x*(1+x*y)*(1-x^2*y)))
sage: plot3d(F,[0,1],[0,1])

EDIT: comments say that as $y\to 0+$, the Hessian becomes indefinite. That is, for my original claim to hold, one need to assume $y>\epsilon>0$, with $\epsilon$ possibly depending upon $x$.

I conjecture a stronger statement: that the function $F(x,y)$ of two variables $x,y$ obtained from $F_n(x)$ by substituting $y$ for $x^{2n}$ is convex for $0<y<x<0$. While this reduces to checking positive definiteness of a matrix of bivariate polynomials of degree about 20, obtained from the Hessian of $F$ by clearing common (positive) denominator, this should be possible to make to work by standard real algebraic geometry tools. Here is an (ugly) plot of $F$:

obtained by exporting the plot from jmol in Sage(math)

sage: var('x y')
sage: F(x,y)=log((x+y^2)*(1+y)*(1-x*y)/(x*(1+x*y)*(1-x^2*y)))
sage: plot3d(F,[0,1],[0,1])

I conjecture a stronger statement: that the function $F(x,y)$ of two variables $x,y$ obtained from $F_n(x)$ by substituting $y$ for $x^{2n}$ is convex for $0<y<x<0$. While this reduces to checking positive definiteness of a matrix of bivariate polynomials of degree about 20, obtained from the Hessian of $F$ by clearing common (positive) denominator, this should be possible to make to work by standard real algebraic geometry tools. Here is an (ugly) plot of $F$:

obtained by exporting the plot from jmol in Sage(math)

sage: var('x y')
sage: F(x,y)=log((x+y^2)*(1+y)*(1-x*y)/(x*(1+x*y)*(1-x^2*y)))
sage: plot3d(F,[0,1],[0,1])

EDIT: comments say that as $y\to 0+$, the Hessian becomes indefinite. That is, for my original claim to hold, one need to assume $y>\epsilon>0$, with $\epsilon$ possibly depending upon $x$.