The best way to solve this is to use Sturm's theorem. This gives an algorithm for computing the number of distinct real roots of any polynomial. The Wikipedia page is quite good, but I'll outline the method here.
Let $f(x)$ be a polynomial. We define a sequence as follows:
$$P_0=f$$
$$P_1=f'$$
$$P_{n+2}=-P_{n}\text{ mod }P_{n+1}$$
where $f'$ is the derivative of the polynomial and, for polynomials $P$ and $Q$, we define $P\text{ mod }Q$ to be the remainder of dividing $P$ by $Q$ - that is, the unique polynomial $R$ of degree less than $\deg Q$ such that $P=cQ+R$ for some other polynomial $c$. (This is also just the result you get by polynomial long division)
For instance, suppose we want to know how many roots $f(x)=x^3+2x+1$ has using this method - of course, we know the answer is $1$, but we should check. We get the following chain:
$$P_0=x^3+2x+1$$
$$P_1=3x^2+2$$
$$P_2=-\frac{4}3x-1$$
$$P_3=\frac{-59}{16}.$$
For any real number $a$, we define $V(a)$ to be the number of sign changes in the sequence $P_0(a),P_1(a),P_2(a),P_3(a)$, where we ignore any zeros. Assuming neither $a$ or $b$ are themselves roots, Sturm's theorem states that $V(a)-V(b)$ is the number of real roots between $a$ and $b$.
Note that $V(-\infty)=\lim_{a\rightarrow-\infty}V(a)$ or $V(\infty)=\lim_{b\rightarrow\infty}V(b)$ are easy to compute by looking at the leading terms of each polynomial. For instance, here we have that $V(-\infty)=2$ since, towards $-\infty$ we have that $P_0$ tends to $-\infty$, $P_1$ to $\infty$, $P_2$ to $\infty$ and $P_3$ is negative - so two sign changes. Then $V(\infty)=1$ because $P_0$ and $P_1$ are positive near $\infty$ and $P_2$ and $P_3$ are negative. This polynomial has $V(-\infty)-V(\infty)=1$ roots, as expected, since it is an increasing function.
This can be a bit laborious to do by hand, but it always works for any polynomial.
The only trick to proving this, at least in the square-free case, is to consider what happens to sign changes in this sequence as one moves along the real line: The number of sign changes can only change near a root of one of the polynomials. However, note that, for some polynomial $c$, we have the following relationship:
$$P_{n}=cP_{n+1}-P_{n+2}$$
Note that if $P_{n+1}$ has a root at a place where $P_n$ doesn't, then near that root, $P_n$ and $P_{n+2}$ must have opposite signs, since $P_n=-P_{n+2}$ at the root. So long as $P_0$ was squarefree (i.e. has no multiple roots), we can note that no consecutive terms share a root, so this always happens. As a result, the zero of $P_{n+1}$ does not affect the number of sign changes. However, if $P_0$ has a root, then the number of sign changes decreases by one there, since, near that root, $f$ and $f'$ have opposite signs prior to the root and equal signs after.
Best Answer
This is a partial answer. I will show that for $a>0$ small enough
$$ f(x) = x(1+x^n)^n - a(1+x^n)^n + a^nx $$
will have exactly one positive real zero and for $a>0$ large enough $f$ will have exactly three positive real zeros.
The case when $a$ is small.
When $a = 0$ the zeros of $f$ are the $n^\text{th}$ roots of $-1$, each with multiplicity $n$, and $x=0$, with multiplicity $1$.
Figure 1: The zeros of $f$ for $n=5$ when $a=0$. The unit circle is shown in blue for reference.
Zeros of polynomials are continuous functions of their coefficients, which means that for any $\epsilon > 0$ and any $n^\text{th}$ root of $-1$, say $\zeta$, $f$ will have exactly $n$ zeros in the disk $|x-\zeta| < \epsilon$ for $a$ sufficiently small. In particular, none of the $n^2$ zeros of $f$ clustered around the $n^\text{th}$ roots of $-1$ will lie on the real axis for $a$ sufficiently small.
Figure 2: The zeros of $f$ for $n=5$ with $0 < a < 0.9$. As $a$ increases from $0$, exactly $5$ zeros erupt from each of the $5^\text{th}$ roots of $-1$. We can also see the zero which was previously located at $x=0$ move to the right. Also, two of the zeros coming from the primitive root of $-1$ and its conjugate travel toward the real axis.
Now $f(0) = -a$ and $f(x) \to +\infty$ as $x \to +\infty$, so $f$ will have at least one positive real zero when $a > 0$. By the above remark this positive real zero will be unique for $a > 0$ sufficiently small.
The case when $a$ is large.
We noted above that $f(0) < 0$ for all $a > 0$. Further, $f(x) \to \infty$ as $a \to \infty$ for fixed $x>0$, which means that $f$ has a positive real zero which tends toward the point $x=0$ as $a \to \infty$.
Numerically we can see that, for $a>0$ large, $f$ has a zero which is approximately equal to $a$ (see Figures 3 and 4). Indeed we see that, for $n>1$,
$$ \begin{align} \frac{f(ax)}{a^{n^2+1}} &= x\left(a^{-n}+x^n\right)^n - \left(a^{-n} + x^n\right)^n + a^{n-n^2}x \\ &\to x^{n^2+1} - x^{n^2} \\ &= x^{n^2}(x-1) \tag{$*$} \end{align} $$
as $a \to \infty$, where the convergence is uniform with respect to $x$ in compact subsets of $\mathbb{C}$. From this we may conclude that $f(ax)$ has a zero arbitrarily close to $1$ for $a$ large enough, and so $f(x)$ has a zero, call it $x_\infty(a)$, for which $x_\infty(a) \sim a$ as $a \to \infty$.
To prove that $x_\infty$ is real we will employ an approach similar to the one we used in Aside 2. Indeed for $a>0$ we have $f(a) = a^{n+1} > 0$ and, for $0 < \epsilon < 1$,
$$ f(\epsilon a) = a (\epsilon - 1) \Bigl(1 + (\epsilon a)^n\Bigr)^n + \epsilon a^{n+1} < 0 \tag{$**$} $$
for $a > 0$ sufficiently large provided $n > 1$. Thus $f$ has a zero in any interval $(\epsilon a,a)$ for $a > 0$ sufficiently large.
Figure 3: The zeros of $f$ for $n=5$ and $1.2 < a < 3$. In addition to the zeros continuing to radiate outwards from the $5^\text{th}$ roots of $-1$ we also see a zero tending toward the point $x=0$ and a large zero quickly moving to the right on the real axis. These are the zeros asymptotic to $a^{1-n}$ and $a$, respectively.
So far we have found two of the $n^2+1$ zeros of $f$. It turns out that the remaining $n^2-1$ zeros of $f$ are asymptotically evenly spaced on the circle of radius $|x| = a^{1/(n+1)}$ as $a \to \infty$.
We have
$$ \begin{align} \frac{f\left(a^{1/(n+1)} x\right)}{a^{n/(n+1)-n-1}} &= a^{-n/(n+1)} x \left(a^{-n/(n+1)} + x^n\right)^n - \left(a^{-n/(n+1)} + x^n\right)^n + x \\ &\to x - x^{n^2} \\ &= x\left(1-x^{n^2-1}\right) \end{align} $$
as $a \to \infty$, where the convergence is uniform with respect to $x$ in compact subsets of $\mathbb{C}$. It follows that $f\left(a^{1/(n+1)} x\right)$ has a zero arbitrarily close to $\omega$ for $a$ large enough, where $\omega$ is any $(n^2-1)^\text{th}$ root of unity. Consequently $f(x)$ has a zero which is asymptotic to $\omega a^{1/(n+1)}$ as $a \to \infty$ for each $\omega$.
In particular, $f$ has a zero which is asymptotic to $a^{1/(n+1)}$ as $a \to \infty$. To show that this is real we proceed as before, observing that $f\left(a^{1/(n+1)}\right) < 0$ and, for any $0 < \epsilon < 1$, $f\left(\epsilon a^{1/(n+1)}\right) > 0$ for $a > 0$ sufficiently large. Thus $f$ has a zero in the interval $\left(\epsilon a^{1/(n+1)}, a^{1/(n+1)}\right)$ for $a>0$ sufficiently large.
Figure 4: The zeros of the rescaled polynomial $f\left(a^{1/(n+1)} x\right)$ for $n=5$ and $0 < a < 5$. We can see the zeros which originally erupted from the $n^\text{th}$ roots of $-1$ (prior to rescaling) travel toward points on the unit circle. This rescaled polynomial also has a large real zero which will tend to $\infty$ as $a \to \infty$.
Figure 5: The last shown state of the zeros displayed in Figure 4. Specifically, these are the zeros of the rescaled polynomial $f\left(a^{1/(n+1)} x\right)$ for $n=5$ and $a = 5$.
From the discussion above we may deduce that the remaining $n^2-1$ zeros do not lie on the real axis for $a>0$ sufficiently large. Thus $f$ has exactly three positive real zeros for $a>0$ sufficiently large.
Summary of the results.
We will suppose always that $a>0$.
For $a$ small the polynomial $f$ has exactly one real zero which is asymptotic to $a$ as $a \to 0$. The remaining $n^2$ zeros tend to the $n^\text{th}$ roots of $-1$ as $a \to 0$.
For $a$ large the polynomial has exactly three real zeros, one asymptotic to $a^{1-n}$, another asymptotic to $a$, and the last asymptotic to $a^{1/(n+1)}$ as $a \to \infty$. The remaining $n^2-1$ zeros are asymptotic to $\omega a^{1/(n+1)},$ where $\omega$ is an $(n^2-1)^\text{th}$ root of unity.