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

@nm  Have you tried a fresh session, and just the simplify commant (no interface(...), restart etc)?

@acer How did you know that

sum( (sum(1/k,k=1..n-1) + 1/n) /n/(n+1), n=1..infinity );

will work? I wouldn't have suspected it!

@nm I meant systematic algorithms. Of course, some empirical algorithms do exist.
Have you tried my example with mma?

@janhardo  A^2 = B^2  does NOT imply A = B.

@sursumCorda PDETools:-Solve works (a bit slower)  without subs too:
PDETools:-Solve(eq,  trial_solution_constants, 'independentof' = t);

@mmcdara Carl has already answered:  (z^(1/n))^n = z by definition of z^(1/n)  for n a nonzero integer.