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(n) &= 1024n^2 - 4096/3n + 320\\
b(n) &= 12288n^3 - 14336n^2 + 4864/3n + 256\\
W(x) &= x^4\left(V(x)-(n+1)^2(2n+1)^2x^{24n-1}(1-x)^4(a(n)x + b(n)(1-x))\right).
\end{align*}

For all $n\ge12$ we claim that $(1-x)^6$ divides $W(x)$, and that all the coefficients of $W(x)/(1-x)^6$ are nonnegative. In particular, $W(x)$ is positive for $x>0$. As $a(n)$ and $b(n)$ is positive for all $n\ge1$, we see that $V(x)$ is positive for all $0<x<1$, provided that $n\ge12$. The cases $2\le n\le 11$ are easily checked directly, for instance by noting that the coefficients of $(1+x)^{\deg V}V(1/(1+x))$ are positive.

I'm not sure about the easiest kind to prove the assertion about the coefficients of $W(x)/(1-x)^6$. We compute them via
\begin{equation}
\frac{W(x)}{(1-x)^6}=W(x)\sum_{k\ge0}\binom{k+5}{5}x^k.
\end{equation}

All what follows is assisted/confirmed by the Sage code below.

The exponents of $W(x)$ have the form $2ni+j$ for $0\le i\le 12$ and $0\le j\le 10$. Let $a_{i,j}$ be the corresponding coefficient. It is a polynomial in $n$.

We have
\begin{equation} 
\frac{W(x)}{(1-x)^6} = \sum_{i,j}a_{i,j}(n)x^{2ni+j}\sum_k\binom{k+5}{5}x^k.
\end{equation}

With $r\ge0$ and $0\le s\le 2n-1$, the coefficient of $x^{2nr+s}$ in $W(x)\frac{1}{(1-x)^6}$
is the sum of all $a_{i,j}(n)\cdot\binom{nr+s-ni-j+y}{5}$ for all pairs $(i,j)$ where $i<r$ and $0\le j\le 10$, or $i=r$ and $j\le s$.

We distinguish the cases $0\le s\le 9$ from the cases $s\ge10$. In these latter cases, we write $s=10+y$. Recall that $s\le 2n-1$, so $10+y=2n-1-o$ for a nonnegative $o$. Upon replacing $n$ with $(11+y+1)/2$, we get polynomials in $n$ and $o$.

In the former cases, the coefficients are polynomials in $n$, which have positive coefficients upon replacing $n$ with $n+12$. So these coefficients are positive for all $n\ge12$. 

In the latter cases it turns out that in all but one case the coefficients are positive. So for nonnegative $y$ and $o$, the values are nonnegative.

For the single exception we see that if we multiply it with a suitable polynomial, the resulting coefficients are positive.

*Remark (answering Jason's question):* It is a known fact that a polynomial $V$ is positive on $(0,1)$ if and only if it is a nonnegative linear combination of polynomials $x^i(1-x)^j$. The problem is that it may involve terms where $i+j$ is bigger than $\deg V$. (This doesn't happen here, though.) I used an LP solver to play a little with $V$ and $\frac{V}{(1-x^2)^4x^{2n-2}}$. From that a pattern showed up which lead to a solution.

By the way, $V$ actually seems to be a nonnegative linear combination of $x^i(1-x)^{\deg V-i}$. This is equivalent to say that all the coefficients of $(1+x)^{\deg V}V(1/(1+x))$ are nonnegative. It is not hard to compute explicit expressions for these coefficients, but I don't see an easy arguments why they can't be negative.

    # 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 = 1024*N^2 - 4096/3*N + 320
    b = 12288*N^3 - 14336*N^2 + 4864/3*N + 256
    W = X^4*(V - (2*N + 1)^2 *(N + 1)^2*X^(24*N-1)*(1-X)^4*(a*X + b*(1-X)))

    # Check that (1-X)^6 divides W(X)
    print all(W.diff(X,i)(X=1).polynomial(QQ) == 0 for i in [0..5])

    # Compute the coefficients a_ij for the exponents 2ni+j of W. 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()]
    K.<n> = QQ[]
    aij = {(i,j):K(0) for i in [0..12] for j in [0..10]}
    for term in l:
        c = term(X=1)                        # get coefficient of term
        e = (term.diff(X)/c)(X=1)            # get exponent of term
        c = K(c.polynomial(QQ))              # convert c to proper polynomial in n
        i, j = ZZ(e.diff(N))//2, ZZ(e(N=0))  # Clumsy method to compute pairs (i,j)
        aij[i,j] += c

    # Check if coefficients aij[i,j] were correctly computed
    Wnew = sum(c(n=N)*X^(2*N*i+j) for (i,j),c in aij.items())
    print (W-Wnew).canonicalize_radical().is_trivial_zero()

    def bino(k): # binomial coefficient binom(-6,k)=binom(k+5,5)
        return prod(k+5-z for z in range(5))/120

    # compute the coefficients of W(X)/(1-X)^6
    K = K.extend_variables(('y', 'o'))
    K.inject_variables()
    for r in [0..12]:
        for s in [0..10]:
            d = 1 if s == 10 else 0
            # compute coefficient of X^(2nr+s+d*y)
            f = sum(aij[i,j]*bino(2*n*(r-i)+s+d*y-j) for i in [0..r]
                    for j in [0..s if i == r else 10])
            # check non-negativity of the polynomial f
            if d == 0:
            # Checks if the coefficient of X^(2nr+s), which is a
            # polynomial in n, has nonnegative coefficients upon replacing
            # n with n+12
                f = f(n=n+12)
                if min(f.coefficients()+[0]) < 0:
                    print "False"
            else:
            # The coefficient of X^(2nr+10+y), which is a polynomial in n
            # and y, has to be nonnegative for 0 <= 10+y <= 2n-1, so 10+y
            # = 2n-1-o for non-negative o. So upon replacing n with
            # (11+y+o)/2, we get a polynomial in y and o which has to be
            # nonnegative for all nonnegative y and o. In all but one
            # case, this holds because the coeffcients are non-negative.
                f = f(n=(11+y+o)/2)
                if min(f.coefficients()+[0]) < 0:
                    c = o^2 + 23*o*y + 1360*y^2 + 99*o + 340*y + 1675
                    print min(c.coefficients()+(c*f).coefficients()) >= 0