If you just want to see the values of i+n as they are computed, you can set printlevel to 2:
restart;
printlevel; #Default value 1
printlevel:=2:
for i from 0 to 10 do
  for n from 0 to 10 do
    i+n
  end do
end do;

##Alternatively, use an explicit print statement, i.e. replace i+n by print(i+n).

restart;
for i from 0 to 10 do
  for n from 0 to 10 do
    print(i+n)
  end do
end do;

I tried
?inert command

where % was mentioned and a link given to value:

?value

I also tried
? %

The page that came up was about the ditto operators, but at the left value (%) was listed as the second entry.

About piecewise.
For the active (non-inert) form to work every other argument has to be a relation or a boolean combination of inequalities, as is stated in the help page for piecewise:
?piecewise

`` is a name, not a relation or of type boolean. The following returns false:

type(``,{relation,boolean});

and this returns true:
type(``,name);

The following version is certainly not simple, but has the advantage of using just one add:

restart;
a:=Matrix([[1,2],[3,4]]):
add(cat('a',b)*cat(f,b),b=indices(a));
evalindets(%,name,eval@parse);
A:=Array(1..2,2..5,0..6,(i,j,k)->i*j*k);
add(cat('A',b)*cat(f,b),b=indices(A));
evalindets(%,name,eval@parse);

I don't know of a way to specify names of arbitrary constants in advance.
But the names can be changed afterwards.
Either by simple assignment e.g.  _C1:=A; _C2:=B; or by the following method which I illustrate in a simple example.

restart;
ode1:=diff(x(t),t)=1+x(t);
ode2:=diff(y(t),t)=2+y(t);
dsolve(ode1); #Uses _C1
dsolve(ode2); #Uses _C1
res:=dsolve({ode1,ode2}); #Uses _C1 and _C2
##Changing names to A and B:
indets(res,name) minus indets({ode1,ode2},name)=~{A,B};
res2:=eval(res,%);

You cannot obtain your desired names by assigning in advance of using dsolve as in _C1:=A; _C2:=B;
because dsolve will not use assigned names.

In the following code I try finding a solution for which phi=0, UF=0, UTF=0.
In view of the difficulty in finding a solution to the original problem, this might be of interest.

I start with your definitions of the equations and the approxsoln GUESS.

systotal:=subs(NBT=NBBT,{ueq2,Teq,eq3,eq4,eq5,eq6,eq7,U(0)=0,UF(0)=0,UT(0)=0,UFT(0)=0,((-cbf*rhobf+cp*rhop)*UF(1)+ rhobf*cbf*U(1))/10000=p2,(((-cbf*rhobf+cp*rhop)*UFT(1)+ rhobf*cbf*UT(1)))/(p2*10000)=T_Bulk,UF(1)=Phiavg*U(1),D(u)(0)=0,u(1)=-lambda*D(u)(1),D(T)(0)=0,phi(0)=phi0,D(T)(1)=1}):
###
sys,bcs:=selectremove(has,systotal,diff):
sys:=fnormal(sys);
#res:=dsolve(sys union bcs,numeric,approxsoln=GUESS,method=bvp[midrich]); #FAILURE

#Setting phi=0, UTF = 0, UF = 0:
eval(sys,{UF=0,UFT=0,phi=0,phi0=0});
sys4:=select(has,%,diff);
eval(bcs,{UF=0,UFT=0,phi=0,phi0=0});
remove(has,%,T_Bulk);
select(hastype,%,function);
bcs4:=remove(has,%,U(1));
nops(bcs4);
GUESS4:=remove(has,GUESS,{T_Bulk,UFT,UF,phi0,phi});
res4:=dsolve(sys4 union bcs4 union {T(0)=0},numeric,approxsoln=GUESS4,method=bvp[midrich]);
plots:-odeplot(res4,[[eta,T(eta)],[eta,u(eta)],[eta,U(eta)],[eta,UT(eta)]],0..1,thickness=3,legend=[T,u,U,UT]);
res4(0);

Solve as a system in one call to dsolve as in

dsolve({eg1,eq2, x(0)=0, ...(fill in inits) },numeric);

But give us your code instead of images so we don't have to type.

Here is another and rather simple minded approach.
It discretizes y and for the integral it simply uses a Riemann integral with fixed interval length d.
It is in need of refinement. But for now here it is:

restart;
eq2:=y(t)=1-h*Int(JacobiTheta3(0,exp(-Pi^2*(t-s)))*y(s)^4,s=0..t); #Better
eq:=y[i]=1-h*Sum(JacobiTheta3(0,exp(-Pi^2*(t[i]-s[j])))*y[j]^4*d,j=1..i-1);
h:=0.65e-4;
d:=.01;
for i from 0 to 100 do t[i]:=i*d; s[i]:=i*d end do:
eqs:=[seq(eq,i=1..100)]:
eqs:=value(%):
res:=fsolve(eqs);
L:=eval([seq([t[i-1],y[i]],i=1..100)],res);
plot(L);

Use an Array instead:

x:=Array(0..2,[1,2,3]);

x[0];

The cause is to be found in the package VectorCalculus. Try the whole worksheet starting with restart and comment out with(VectorCalculus) and BasisFormat(false).

Then the double summation gives you 3.

I use this package very rarely and, when I do, I use the long form as in VectorCalculus:-Jacobian.

Notice that the package redefines `*`,`+`,`-`,`.`,<,>,<|>  and many more familiar operations or procedures.

I found one solution (so far):

sol:=dsolve(newsys2 union bcs union {bcs2,(D@@2)(h2)(1)=1}, numeric,approxsoln = [omega2 = 2.08*10^13, h1(theta) = 5*10^7*theta*(1-theta)*(1/2-theta), h2(theta) = 0.05, h3(theta) = 7*10^6*theta*(1-theta)*(1/2-theta), h4(theta) = 5*10^6*theta*(1-theta)*(1/2-theta)], abserr = 0.1,initmesh=1024);

plots:-odeplot(sol,[seq([theta,cat(h,i)(theta)],i=1..4)],thickness=3);



h3 is much smaller than the other three.
plots:-odeplot(sol,[theta,h2(theta)]);


The corresponding omega2 is found from sol(0) to be omega2 = 2.0798272963267560997*10^13. omega is the square root of that.

To find a successful approxsoln I used my own homemade procedure.
Since the system is linear and homogeneous any constant multiple k of (h1,h2,h3,h4) will also be a solution to the ode system for the same omega2. However, the inhomogeneous boundary condition (D@@2)(h2)(1) = 1 wiil have to be adjusted to (D@@2)(h2)(1) = k.
And indeed the following, where the factor k is 10^(-7) also works:
dsolve(newsys2 union bcs union {bcs2,(D@@2)(h2)(1)=10^(-7)}, numeric,approxsoln = [omega2 = 2.08*10^13, h1(theta) = 5*theta*(1-theta)*(1/2-theta), h2(theta) = 0, h3(theta) = 0.7*theta*(1-theta)*(1/2-theta), h4(theta) = 0.5*theta*(1-theta)*(1/2-theta)], abserr = 0.1,initmesh=1024);

You are missing a few multiplication signs.
There are 3 solutions.

restart;
eq1:=y=86:
eq2:=y=-0.0000054527*x^3+0.010903836*x^2+0.0714244709*x+74.18816:
sol:=solve({eq1,eq2},{x,y});
##Graphical check:
plot([rhs(eq1),rhs(eq2)],x=-37..31);
plot([rhs(eq1),rhs(eq2)],x=-37..2010);
plot([rhs(eq1),rhs(eq2)],x=-37..2010,y=0..90);
plot(ln~([rhs(eq1),rhs(eq2)]),x=-37..2010); #The logarithm applied to both

You could do:

plot([[t^2 + t + 41,t^2 + 40,t=-15..15],[4*t^2 + 163,2*t^2 + t + 81,t=-7..7]]);

I'm assuming that when you differentiate w.r.t. t you actually mean x?
There is no parameter in the system, so I don't quite understand the fitting part.
Whatever you get out of solving the system is it.
There is an immediate problem with T(0)=0 since T(x) appears in the denominator. That problem is probably minor since the problem is in exp(-10/T(x)), which has limit 0 when T(x)->0 from the right.
Your system can be reduced to just an ode in one of Y or T. I chose T.

If you insist on the initial conditions as they are, then clearly
T(x) = 0, Y(x) = 1 for all x is the solution (using that exp will be zero as mentioned).

So you ought to rethink the problem.

Finally, I doubt that you will get a symbolic answer. Certainly it doesn't come out of dsolve right away.

restart;
sys:=diff(Y(x),x,x)=500000*Y(x)*exp(-10/T(x))+diff(Y(x),x),diff(T(x),x,x)=10*(-100000*Y(x)*exp(-10/T(x)))+diff(T(x),x);
ode:=2*sys[1]+sys[2];
icsY:=Y(0)=1 , D(Y)(0)=0; icsT:= T(0)=0 , D(T)(0)=0;
dsolve({ode,icsY},Y(x));
Yres:=eval(%,{icsT});
odeT:=subs(Yres,sys[2]);
eval(odeT,{exp(-10/T(x))=0,T(x)=0});
eval(Yres,{T(x)=0,Y(x)=1});

You clearly have too many boundary conditions, so I interpret your problem as one of determining omega so that the extra boundary condition is satisfied.
You will have to use numerical solution.
(Pi should be spelled like that, NOT pi).

restart;
ode:=diff(y(x), x, x) = -(8*omega*(1-exp(-8*x))*exp(-8*x)/(2*x^2)-Pi^2/(32*x))*y(x);
bcs:=y(0)=0,D(y)(0)=1,y(1/2)=0;
dsolve(ode); #We see that no symbolic solution is to be expected
res:=dsolve({ode,bcs},numeric,method=bvp[middefer]);
plots:-odeplot(res,[x,y(x)],0..1/2);
res(0); # You see that omega=0.805476451353708

Three various ways using elementwise operations (~) .

L:=[-3, -2.5, 0, 1, 3/2, Pi, sqrt(5)];
type~(L,nonnegint); #Just the results in order
L=~type~(L,nonnegint); #Equations
zip(`[]`,L,type~(L,nonnegint)); #Lists