restart: assume(theta, integer): Digits := 100: 
with(RegularChains): with(SemiAlgebraicSetTools): with(FastArithmeticTools): 
with(Groebner);interface('prettyprint' = 0): 

c := proc (e) options operator, arrow: e^theta end proc:       ### cost function
f := proc (c) options operator, arrow: c^(1/theta) end proc:  ### inverse of cost function

U1 := proc (k) options operator, arrow: sum(q[j]*(v1-c(f(v2)+y[j]-x[k]-d)), j = 1 .. k)-c(x[k]) end proc: 
dU1 := proc (k) options operator, arrow: simplify(diff(U1(k), x[k])) end proc: 
U2 := proc (k) options operator, arrow: sum(p[j]*(v2-c(f(v1)+x[j]-y[k]+d)), j = 1 .. k-1)-c(y[k]) end proc: 
dU2 := proc (k) options operator, arrow: simplify(diff(U2(k), y[k])) end proc: 
eqU1 := proc (n) options operator, arrow: seq(U1(k)-u[1], k = 1 .. n) end proc: 
eqU2 := proc (n) options operator, arrow: seq(U2(k), k = 2 .. n+1) end proc: 
eqdU1 := proc (n) options operator, arrow: seq(dU1(k), k = 1 .. n) end proc: 
eqdU2 := proc (n) options operator, arrow: seq(dU2(k), k = 2 .. n+1) end proc: 
prob := proc (n) options operator, arrow: {1-(sum(p[j], j = 1 .. n)), 1-(sum(q[j], j = 1 .. n+1))} end proc: 

C := proc (a) options operator, arrow: subsindets(a, indexed, convert, symbol) end proc:
CC:=proc(a): subsindets(a,symbol,parse);end proc:

aux := proc (n) options operator, arrow: [eqU1(n), eqdU1(n), eqU2(n), eqdU2(n), op(prob(n))] end proc:
sys := proc (n, t) C(subs(d = 0, v1 = 1, v2 = 1, theta = t, y[1] = 0, aux(n))) end proc:
v := proc (n) local z, i: z := []: for i to n do z := [y[i+1], x[i], op(z)] end do: z; end proc: 
vars:= proc (n) local z, i: z := []: for i to n-2 do z := [p[i], q[i], op(z)] end do: 
op(C([p[n], q[n], p[n-1], u[1], q[n-1], op(z),op(v(n))])) end proc: 

RR := proc (n): PolynomialRing([vars(n)]) end proc:

PP := proc (n) local i, z, w: z := v(n): for i from 1 to nops(z)-1 do z[i] := z[i]-z[i+1] end do: 
z := z[1 .. nops(z)-1]: C([seq(p[i],i=1..n),seq(q[i],i=1..n+1),1-u[1],1-y[n+1], op(z),x[1]]) end proc:

i:=5;
F0:= sys(i-1,2);
F1:= sys(i,2):
v0:=SuggestVariableOrder(F0);
v1:=SuggestVariableOrder(F1);
N := [C(u[1])]: H := []: P1 := PP(i):
R1 := PolynomialRing([v1]):
B0:=Basis(F0,plex(v0)):
B1:=Basis(F1,plex(v1)):
Q0:=simplify(B1[1]/B0[1]);
Q1:=[Q0,op(B1[2..nops(B1)])]:
S:=RealRootIsolate(Q1,N,P1,H,R1,method='Discoverer');
Display(S,R1);
quit;