When you do

dsolve({Teq,UEQ});

you get a result that gives you 3 different solutions for T(eta):

1. T(eta) = 0, which you cannot use because of your boundary conditions
2. T(eta) = _C3*eta+_C4
3. A complicated looking expression giving T(eta) in terms of a RootOf. Out of this one you are surely not going to get an analytic answer.

u(eta) is given as an unevaluated double integral involving T(eta).

So the best (or only) hope is for T(eta) = _C3*eta+_C4 to give us an expression for u(eta) where the integrals can be computed.
By using the boundary conditions we get T(eta) = eta.
After that we try

dsolve(eval(UEQ,T(eta)=eta));

We notice that we no longer have a double integral, but an integral, which Maple cannot do as it stands:
value(%);

You can try
dsolve({eval(UEQ,T(eta)=eta),D(u)(1)=0,u(0)=L*D(u)(0)}); #You don't get anything
dsolve({eval(UEQ,T(eta)=eta),u(0)=L*D(u)(0)}); #You do get something but an integral is still there
dsolve({eval(UEQ,T(eta)=eta),D(u)(1)=0}); #Also here you get an integral
############
## Numerical solution:
gama1:=0.2:
rhop:=5180:
rhobf:=998.2:
a[mu1]:=39.11:
b[mu1]:=533.9:
a[k1]:=7.47:
b[k1]:=0:
L:=1: N_bt:=1: #Not supplied by you.
res:=dsolve({Teq,UEQ,D(u)(1)=0,u(0)=L*D(u)(0),T(1)=1,T(0)=0},numeric);
plots:-odeplot(res,[eta,T(eta)-eta],0..1);
## In this case we do get a solution where T(eta) = eta. See note below ******
#This gives us not much of a hint as to the form of the value of the integral for u:
plots:-odeplot(res,[eta,u(eta)],0..1);
res(0);
res(1);
####
***Note: This is not surprising in view of the result of
`dsolve/numeric/bvp/convertsys`({Teq,UEQ},{D(u)(1)=0,u(0)=L*D(u)(0),T(1)=1,T(0)=0},[T,u],eta,l);
Or the result of:
solve(Teq,diff(T(eta),eta,eta));

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