Algebraic Question

The following system comes up when analyzing one generalization of the Rosenbrock function.
	0==(Subscript[x, i]) (-2 (Subscript[x, i+1])+ϵ+1)+2 (Subscript[x, i]^3)-(Subscript[x, i-1]^2)-ϵ
Numerical evidence suggests that there is omse interesting structure to the real solutions of these "circulant" polynomials.  Circulant because I want to identify the n+1 st variable with the 1st variable for some value of n and the 0th varaible with the nth.  The range of ϵ values that are of interest are 0≤ϵ<1/8.  The particular values of relevance for the most common Rosenbrock functions are 1/100, 1/105 and 1/104 so we are particularly interested in ϵ  fairly small. The solution structure does not seem to depend all that much on the value of ϵ:  I have no idea why!

1) There are three distinct real solutions for n=2,3.  They are constant i.e. (Subscript[x, i])==(Subscript[x, i+1]).
One of the solutions is the trivial {1,1,…,1}.  The other is near 0.  The third near 1/2.
These values can be expressed analytically  {1,(1/4)(1±(Sqrt[1-8ϵ])}.
2) This pattern seems to continue for n=4,5,6
3) For n=9 I thought there was a curve of solutions but I now believe that it is two families of solutions linked by the shift invariance. See the images below. 
4) For n=7,8  there does not seem to be a similar curve of solutions.  I do not think Mathematica's analytical algorithm will finish for n=8 in any reasonable amount of time.
5) For n=8 and up I can use a numerical solver to find accurate numerical solutions.

n=4;
Eqns=Table[
0==(Subscript[x, i]) (-2 (Subscript[x, i+1])+ϵ+1)+2 (Subscript[x, i]^3)-(Subscript[x, i-1]^2)-ϵ,
{i,1,n}]/.{(Subscript[x, 0])->(Subscript[x, n]),(Subscript[x, n+1])->(Subscript[x, 1])}
Vars=Table[(Subscript[x, i]),{i,1,n}]
Tol=(10^(-6));
Select[
Vars/.NSolve[Eqns/.ϵ->0.01,Vars],
Max[Abs[Im[#]]]<Tol&
]

The numerical solver looks like

n=9;
Eqns=Table[
0==(Subscript[x, i]) (-2 (Subscript[x, i+1])+ϵ+1)+2 (Subscript[x, i]^3)-(Subscript[x, i-1]^2)-ϵ,
{i,1,n}]/.{(Subscript[x, 0])->(Subscript[x, n]),(Subscript[x, n+1])->(Subscript[x, 1])}
Vars=Table[(Subscript[x, i]),{i,1,n}]

NumSol=134;
Data=Table[
Vars/.FindRoot[
Eqns/.ϵ->0.01,
Transpose[{Vars,RandomReal[{0.2,0.8},n]}],
MaxIterations->(10^4),
PrecisionGoal->9
],
{NumSol}];

There seems to be a solutions with lots of different intial (Subscript[x, 1])values.

SortedData=Sort[Data];
x1Vals=SortedData[[All,1]];
ListPlot[x1Vals,PlotRange->All]

To visualize the solution curve I will plot these guys as ordered pairs {(Subscript[x, i-1]),(Subscript[x, i])}

ListLinePlot[
NineData=Table[
Temp=Data[[i]];
Transpose[ {Temp,RotateRight[Temp]}],
{i,Length[Data]}],
PlotLabel->{n,0.01},
PlotRange->All
]

I believe the inside ones are good solutions.  I believe the outside ones are not! but I could be wrong!  I believe, the inside oval is actually two distinct solution sets linked by the shift invariance.

ListPlot[
Take[NineData,14],
PlotLabel->{n,0.01},
PlotRange->All
]

Here is the question.

"What, if anything, can be said about the number of real distinct solutions to such circulant polynomial systems?"