@Andiguys Your C_1 is
{0 <= 0.000023339978559*p1 - 0.85357635872}

Probably you want to replace the previous plot2 with:

plot2 := plot([solve(rhs(C_1[])), t, t = -0.2e9 .. 0.1e9])

I don't understand why you complicate things for such a simple problem.

@Andiguys 

plot1:=plot(TRC(p1), p1 = 0 .. 60000):
plot2:=plot(rhs(C_1[]), p1=10000..60000, thickness=3,color=blue):
plots:-display(plot1,plot2);

Not very useful because C_1 has a small slope.

But you don't need it as I said.

@Andiguys You have a simple concave function TRC in one variable (p1);  see  simplify(TRC(p1)).
You can compute the exact max using

maximize(TRC(p1), location);

                  33379441.96, {[{p1 = 24171.82839}, 33379441.96]}

If you want the max value in another interval [a,b], it will be   33379441.96 if  24171.82839 is in [a,b]
and max(TRC(a), TRC(b)) otherwise.

@Carl Love  Unfortunately, Maple is not able to compute inf, liminf, sup,limsup [Mathematica does!]. When the limit of a bounded function does not exist, limit could return a range a..b including the the image of the function.

For example,

limit(sin(n)+cos(n), n=infinity);

returns -2..2  but liminf and limsup are -sqrt(2), sqrt(2)  (even for sequences).

@Preben Alsholm Yes, as it should, but the fact that assuming makes it fail must be investigated.

@acer  simplify fails if  assuming n::posint   is appended!!!

@Carl Love The simplified result is (6*n - 5) / 16.
It is frustrating that Maple needs a lot of help to obtain it.

@Preben Alsholm  The maplesoft download page (www.maplesoft.com/support/downloads/) contains the update.

Does the interrupt button (bottom left) work in 2025.1? Try:

a:=nextprime(10^83):
b:=nextprime(10^95):
x:=a*b:
ifactor(x);

@mmcdara  I had the impression that you are interested in the sum of (random) floats.

For the symbolic case, see: evalf random results - serious bug - MaplePrimes

@C_R   The main problem is that if a() is defined from the beginning, Maple gives the wrong solution (Carl's code is a workaround).

restart;
a := z-> 1+2*Heaviside(z):
de:= Diff(a(x)*Diff(u(x),x),x) = -1:
bc:= u(-1)=0, u(1)=0: 
sol:= dsolve({de,bc}, u(x));

@acer I see that many users here are interested in a rather low level maths but with nice typesetting.
For such purpose, a better software used to be Scientific Workplace; it produced a good LaTeX output.  Actually, it contained a Maple engine as a CAS, which was later replaced with MuPAD. Unfortunately, Scientific Workplace is now discontinued.

@Rouben Rostamian  The problem is that dsolve works AFAIK at formal level and does not really care about generalized (weak) solutions.

You could even rise your expectations:

A:=(((sin(sqrt(3)/2)*sqrt(3) - 3*cos(sqrt(3)/2))*cos(sqrt(3)*x/2) - 
sin(sqrt(3)*x/2)*(sqrt(3)*cos(sqrt(3)/2) + 3*sin(sqrt(3)/2)))*exp(-1/2 + x/2))/3 :
combine(A):
convert(%, phaseamp, x);

@dharr Very efficient. Vote up!