The following (hopefully) settles the question. Fix $n$, and let $V(x)$ be the polynomial from the question, for which we want to show that it is positive on the open interval $(0,1)$. Set \begin{align*} a &= 64n^2 - 256/3n + 20\\ b &= 768n^3 - 1024n^2 + 272n - 24\\ c &= 14336/3n^4 - 306304/45n^3 + 35792/15n^2 - 6104/15n + 32 \end{align*} Then \begin{equation} \frac{V(x)}{(1-x^2)^4x^{2n-2}}=(1-x)^3H(x)+(n+1)^2(2n+1)^2x^{22n-4} (ax^2+bx(1-x)+c(1-x)^2) \end{equation} for a polynomial $H(x)$.
Here $a,b,c$ are positive for $n\ge2$. We show that all coefficients of $H(x)$ are positive whenever $n\ge6$.
I'm not sure about the easiest kind to prove the assertion about the coefficients of $H(x)$. But the following works: Multiplying the equation by $(1-x)^4(1+x)^7$ yields \begin{equation} (1-x^2)^7H(x)=W(x) \end{equation} where $W(x)$ is a polynomial for which we have an explicit (though messy) expression. Write $H(x)=H_0(x^2)+xH_1(x^2)$ and $W(x)=W_0(x^2)+xW_1(x^2)$. Then $H_0(x)(1-x)^7=W_0(x)$ and $H_1(x)(1-x)^7=W_1(x)$. Thus we are reduced to show that $H_0(x)$ and $H_1(x)$ have nonnegative coefficients. We obtain these coefficients by multiplying $W_\star(x)$ with the formal power series $\frac{1}{(1-x)^7}=\sum\binom{k+6}{6}x^k$.
All what follows is assisted/confirmed by the Sage code below.
There is a finite set of pairs $(i,j)$ of nonnegative integers, independent of $n$, such that the exponents of $W(x)$ have the form $2ni+j$. Let $a_{i,j}$ be the corresponding coefficient. It is a polynomial in $n$. It turns out that $0\le i\le 11$ and $-4\le j\le 11$.
We have \begin{equation} W = \sum_{i,j}a_{i,j}(n)x^{2ni+j}. \end{equation} Thus we obtain $W_0$ or $W_1$ according as to whether $j$ is even or odd.
The exponents of $W_0(x)$ now are $ni+j/2$ for even $j$, and those of $W_1(x)$ are $ni+(j-1)/2$ for odd $j$. For better readability, we call the modified $j$'s again $j$. So the exponents are $ni+j$ for certain pairs $(i,j)$, and the associated coefficients will be called $a_{i,j}$. After this modification, we have $-2\le j\le 5$.
If $n\ge8$, then $ni+j>ni'+j'$ if and only if $i>i'$ or $i=i'$ and $j>j'$. This gives a total order on the pairs $(i,j)$. Every exponent of $x$ in $W_\star(x)/(1-x)^7$ which matters has the form $ni+j+y$ where $0\le y<n(i'-i)+j'-j$, where $(i',j')$ is the successor of $(i,j)$ for our order.
The coefficient of $x^{ni_1+j_1+y}$ in $W_\star(x)\frac{1}{(1-x)^7}$ is the sum of all $a_{i_1,j_1}(n)\cdot\binom{ni+j-ni_1-j_1+y+6}{6}$ for all $(i_1,j_1)\le(i,j)$.
Thus we obtain finitely many (actually $126$) explicit polynomials $(n,y)$ in $n$ and $y$ (of small degrees in $n$ and $y$) for which we have to prove that they are positive for all $n\ge8$ and all $y$ in the specified range.
We distinguish two cases: If $(n'i+j')-(ni+j)=1$ ($(i',j')$ the successor of $(i,j)$ again), then only $y=0$ is possible. In all cases, we see that $f(n+6,0)$ has nonnegative coefficients. Thus $f(n,0)\ge0$ whenever $n\ge6$.
Now assume that $(n'i+j')-(ni+j)>1$. In these cases it turns out that $i'=i+1$, so $0\le y< n-e$ for some integer. Thus $y=n-e-1-o$ for some non-negative $o$. Thus $f(n,o)=f(y+e+1-o,o)=g(y,o)$. This latter polynomial, in $y$ and $o$, has to be non-negative for all non-negative $o$ and $y$. It turns out that except for a single example actually oll the coefficients are nonnegative. For this single exception, $g(y,o+1)$ has nonnegative coefficients. Thus $g(y,o)\ge0$ except possibly for $o=0$. But $g(y,0)$ has nonnegative coefficients, and things are settled.
Our argument used $n\ge8$. For $n=6$ and $7$, one can directly compute $H(x)$ and check that the coefficients are nonnegative. The cases $2\le n\le5$ are checked directly for $V(x)$.
Here is the Sage code:
# Formally compute W
var('x n')
P = (1+x^(4*n-1))*(1+x^(2*n))*(1-x^(2*n+1))
Q = (1+x^(2*n+1))*(1-x^(2*n+2))
V = P*Q^2*P.diff(x,2)+(P*Q.diff(x))^2-P^2*Q*Q.diff(x,2)-(P.diff(x)*Q)^2
a = 64*n^2 - 256/3*n + 20
b = 768*n^3 - 1024*n^2 + 272*n - 24
c = 14336/3*n^4 - 306304/45*n^3 + 35792/15*n^2 - 6104/15*n + 32
W = V*(1+x)^3*x^(2-2*n)-(1-x)^4*(1+x)^7*(n+1)^2*(2*n+1)^2*x^(22*n-4)*(a*x^2 + b*x*(1-x)+ c*(1-x)^2)
# Compute the pairs (i,j) for the exponents 2ni+j of W, and the
# corresponding coefficients a_ij. Somewhat clumsy, as I don't know how
# to deal with polynomials where exponents are symbolic expressions
#
# l is the list of summands of W
l = [z.canonicalize_radical() for z in W.expand().operands()]
aij = {}
for term in l:
c = term(x=1) # get coefficient of term
e = (term.diff(x)/c)(x=1) # get exponent of term
i,j = ZZ(e.diff(n))//2,ZZ(e(n=0)) # Clumsy method to compute pairs (i,j)
aij[i,j] = aij[i,j] + c if (i,j) in aij else c
# Check if the pairs (i,j) and coefficients aij[i,j] were correctly computed
Wnew = sum(c*x^(2*n*i+j) for (i,j),c in aij.items())
print (W-Wnew).canonicalize_radical().is_trivial_zero()
# pairs (i,j) and coeffcients for W_0 and W_1, respectively
aij0 = {(z[0],z[1]//2):c for z,c in aij.items() if is_even(z[1])}
aij1 = {(z[0],(z[1]-1)//2):c for z,c in aij.items() if is_odd(z[1])}
# compute and collect coefficients of W_0(x)/(1-x)^7 and W_1(x)/(1-x)^7
var('y o')
hcoeffs = []
for aij in [aij0,aij1]:
aijk = aij.keys()
aijk.sort()
for k in range(len(aijk)-1):
(i,j) = aijk[k]
# f(n,y) is the coefficient of x^(ni+j+y) for 0 = 1.
print min(f(o=o+1).polynomial(QQ).coefficients()) >= 0
# So only o=0 remains to check. We see that f(o=0) has nonnegative coeffcients
print min(f(o=0).coefficients()) >= 0