f:=unapply(simplify(sum(i, i=1..n)), n);
f(n);
f(10);

Unfortunately there are no solutions to this system that completely satisfy your conditions. First we find all the solutions of this system (real and imaginary). Maple easily solves polynomial systems, so first we squared the last 2 equations. The resulting system will be equivalent to the original system under conditions  x>=0, x+y>=0 .

restart;
eq1:=4.1*10^9-x*Tl^4-4.527587408*10^7/Pi=0:
eq2:=4.1*10^9-(0.6013906720-x-y)*Th^4+4.527587408*10^7/Pi=0:
eq3:=4.1*10^9-Ti^4*y=0:
eq4:=1-1.487209994*10^(-20)*Tl^8/(1+sin(-0.41)^2)=x^2:
eq5:=1-1.487209994*10^(-20)*Th^8/(1+sin(-0.41)^2)=(x+y)^2:
_MaxSols:=infinity:
sol:={solve({eq1,eq2,eq3,eq4,eq5}, {x,y,Tl,Th,Ti})}:
sol1:=simplify(sol, zero):
nops(sol1);

                                                              1024

Maple found 1024 solutions, among which there are many imaginary solutions. Next we first try to take all the real solutions and see that there are no such solutions where all the roots are real. But there are solutions where all the roots are real except for Ti :

S1:=select(t->eval(`and`(Im(x)=0, Im(y)=0, Im(Tl)=0, Im(Th)=0, Im(Ti)=0), t), sol1);
S2:=select(t->eval(`and`(Im(x)=0, Im(y)=0, Im(Tl)=0, Im(Th)=0), t), sol1):
nops(S2);
S3:=select(t->eval(`and`(x>=0, x+y>=0, Tl>Th), t), S2);
nops(S3);

This can be done in a mass of ways. If you want the cone to have the specified proportions, you first need to do some calculations to find out the necessary parameters (vertex coordinates, base radius and height).

The first (easy) way:

A:=plottools:-cone([0,0,-1], 1, 1, style=surface, color="Green"):
B:=plottools:-sphere(color=pink, transparency=0.5):
C:=plottools:-reflect(A, [[0,0,0],[1,0,0],[0,1,0]]):
plots:-display(C, B, axes=normal, view=[-1.4..1.4,-1.4..1.4, -1.2..1.4], scaling=constrained);

Addition. The second way.

 Let us  inscribe a cone , whose height is twice the radius of the base, in the unit sphere. First we need to find the radius of the base of this cone:

R:=solve(2*x=sqrt(1-x^2)+1);  # Radius of the base of the cone
A:=plot3d([r*cos(t), r*sin(t), -sqrt(1-R^2)], r=0..R, t=0..2*Pi, style=surface, color="Green"):  # The base of the cone
B:=plot3d([r*cos(t), r*sin(t), -2*r+1], r=0..R, t=0..2*Pi, style=surface, color="Green"): # Lateral surface of the cone
C:=plottools:-sphere(color=pink, transparency=0.7):
plots:-display(A, B, C, axes=normal, view=[-1.2..1.2,-1.2..1.2, -1.2..1.4], scaling=constrained);

Here the cone surfaces are given in terms of their parametric equations.

I supposed that  A, B, alpha, y  are reals,  i=I  is imaginary unit.

restart;
w := A*exp(-alfa*(1+I)*y)+B*exp((1+I)*y);            
u := evalc(w);
subs(I=0, u);

eq:=diff(y(x),x)=sqrt(1+a*x+2*y(x)):
eq1:=eval(eq, a=0);  # Substitution a=0 into the equation eq
sol:=dsolve({eq1, y(0)=1}, y(x));  # Solution of the equation eq1
eval(eq1, sol);  # Substitution of the solution into the equation eq1
is(%) assuming x>0;  # Check of the solution
eval(sol, x=1);  # Exact value y(1)
evalf[5](%);  # Approximate value y(1)

F:=unapply(int(1/(0.1060444429e-1-0.2120888857e-1*X+0.1033933318e-1*X^2), X), X);
fsolve(F(x)-F(0)=5000, x);
simplify(fnormal(%), zero);

See  LinearAlgebra:-GenerateMatrix  command in Help.

Expr:=-(1/4)*G*r^2/eta-(1/4)*G*(ri^2-ro^2)*ln(r)/(eta*(ln(ri)-ln(ro)))-(1/4)*G*(ln(ri)*ro^2-ri^2*ln(ro))/(eta*(ln(ri)-ln(ro))):
Fact:=-(1/4)*G/eta:
``(Fact)*map(t->t/Fact, Expr); 

Addition.  This method does not work if the initial expression has members that do not contain the common factor to be taken out (below I have marked these members through C). Here's a way of solving the problem:

restart;
Expr:=-(1/4)*G*r^2/eta-(1/4)*G*(ri^2-ro^2)*ln(r)/(eta*(ln(ri)-ln(ro)))-(1/4)*G*(ln(ri)*ro^2-ri^2*ln(ro))/(eta*(ln(ri)-ln(ro)))+C:
Fact:=-(1/4)*G/eta:
solve(t=Fact, G):
subs(t=``(Fact), collect(subs(G=%, Expr), t));           

as an example:

restart;
sol := dsolve({diff(i(t),t)+i(t)=0, i(0)=2}, i(t));
assign(sol);
i:=unapply(i(t), t);
plot(i(t), t = -1 .. 4, view = -1 .. 3, gridlines);
i(1);

Edit.

A := day(hello(a, b), c): 
op(0, A);
op([1, 0], A);
                                                                day
                                                               hello

Addition. If you want to extract the names of all functions simultaneously, you can use the following procedure:

restart;
P:=proc(Expr)
if op(Expr)=Expr then return NULL else
op(0,Expr), map(t->P(t), [op(1..-1, Expr)])[ ] fi;
end proc:

Example of use:

A := day(hello(a, morning(b)), night(c)):
P(A);
                                           day, hello, morning, night

Edit.

restart;
a:=convert(-1.327553040*10^37, fraction);
# Obviously the given number is 10^b where
b:=evalf[50](a/ln(10));
n:=floor(b);
c:=evalf[50](b-n);
10^c*10^``(n);   # The answer in scientific notation with standard 10 digits accuracy

First, you can solve the equation for each value  T  in the range  300..800  with step 1, and then plot a graph using the data obtained (you will get a nice plot of even better quality than lower):

restart;
Data:=[seq([t,fsolve(eval(x*(3-2*x)^2/4/(1-x)^3=18325.73901+exp(36.58420421-10902.09286/T), T=t))], t=300..800)];
plot(Data);
 

Before the plotting both parts of the equation are multiplied by (1-x)^3 :

restart;
plots:-implicitplot(x*(3-2*x)^2/4=18325.73901*(1-x)^3+exp(36.58420421-10902.09286/T)*(1-x)^3, T=300..800, x=0..1, gridrefine=5);

In fact, this line intersects the circle at two points. Here are 2 ways to isolate the coordinates of these points (the second way is similar tomleslie's one only for two solutions):

eq1:= x^2+y^2=0.314:
eq2:= y=0.05180967688*x:
Sol:=solve({eq1,eq2}, [x,y]);
P1:=[seq(eval([x,y], s), s=Sol)];   # The first way
P2:=map(rhs~, Sol);   # The second way

In Help see  DocumentTools:-Tabulate  command.