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$:

[![plot on [0,1]x[0,1]][1]][1] 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$.  

  [1]: https://i.sstatic.net/l9j8a.jpg