@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!

@Rouben Rostamian  Nice analysis, vote up!

@Andiguys The variable a is already in use. To keep it alive, rename the new one e.g.
aa := lhs(C_1[])
In NLPSolve too.

@Andiguys 
As my example shows, you must use:

a:=lhs(C_1[]); 
b:=rhs(C_2[]);
NLPSolve(TRC(Pc), Pc = a..b, maximize, method=branchandbound);

                        a := 256.9131467
                        b := 1436.173707
             [9.87515033041729*10^8 , [Pc = 1436.17370700000]]

@Andiguys In my answer, P is the numerator of C and is simplified (actually normalized) simbolically.

@janhardo No, I have obtained a very similar answer via Copilot with "Think Deeper"" option.

@Alfred_F  You should be aware that the "Maple solution" is just typesetting. The only computation made by Maple is Zeta(2)^2/Zeta(4) = 5/2.  The pictures at the end are screenshots generated by AI, probably Microsoft Copilot.

@Alfred_F The standard method is to use the Moebius inversion. More generally, summing 1/(a^p*b^p), over (a,b)=1, the result is Zeta(p)^2 / Zeta(2*p).

@mmcdara Vote up of course!
I have checked some Maple versions >= 2018 I have installed, but they all do not work!
It would be nice to know when and what has changed (in sum, limit) in the interval 2015--2018.

@Alfred_F For such unusual problems you must guide Maple as I did.
BTW, Mathematica is also unable to compute the limit directly.