Download 2curves-vv.mw
(edited for Global max).

The constraints must be included in a set.

Optimization:-NLPSolve( TP_T(p1, p2), {C1, C2}, assume = nonnegative,  maximize );

diff( c * x / (1+x), x);

Optimization:-Maxima computes only local maxima.
In your case when computing the max in the interval {256.9131467 <= Pc, Pc <= 1436.173707},
Pc = 256.9131467  is a local maximal point. The global maximal point is of course Pc = 1436.173707.
If you want a global max (1-dimensiomal only!) you must use NLPSolve (and also method=branchandbound,   but it seems to be the default here).

When you meet such an aparent discrepancy you should look at a simple example; it is also much better to post a minimal example and with standard notations. Here is such an example:

f:=(x-2)^2:
Optimization:-Maximize(f, {1 <= x, x <= 5});
Optimization:-NLPSolve(f, x=1..5, maximize); # method=branchandbound);
plot(f, x=0..5);

                         [1., [x = 1.]]

                         [9., [x = 5.]]

C := C1 - C2 is a rational function (with 15 indeterminates). It will be enough to compute the minimal value for C and - C.
This can be obtained with the command  DirectSearch:-GlobalOptima(C, vars);
Here DirectSearch is a package which can be downloaded at https://www.maplesoft.com/Applications/Detail.aspx?id=101333

Actually, both min values are < 0, so, the answer to your question is that C1<C2 and C2<C1  are both false.
This can be seen directly (without DirectSearch) using:

P:=numer(normal(C1-C2)) / (lambda*varphi):
eval(-P, {Cm = 1, Cr = 10, Pu = 1, U = 2, a = 9, beta = 9, k = 10, lambda = 10, p1 = 1, tau0 = 10, upsilon = 10, varphi = 10, w = 10, varepsilon = 10, U[0] = 10});
#                     -63035971982554554000

eval(P, {Cm = 89, Cr = 65, Pu = 81, U = 50, a = 37, beta = 43, k = 44, lambda = 0, p1 = 85, tau0 = 74, upsilon = 75, varphi = 59, w = 62, varepsilon = 48, U[0] = 41});
#                       -19372211411803488

Download LM-vv.mw

Maple cannot compute the limit, but your expression 
Sum(Sum(j/(j^2 + k^2), j = 1..n), k = 1..n) / n
is a Riemann sum of the double integral

int(y/(x^2+y^2), [x=0..1,y=0..1]);

and this is the result of the limit (not difficult to justify).

Download test12-vv.mw

B:=(4*exp(x/2) - 3*x - 2*exp(x/3)+1)^11 + (4*exp(x/2) + x - 2*exp(x/3)+1)^11 - 13:
A:=combine(expand(B)):
A1:=simplify(A, size): B1:=simplify(B, size):
length~([A, B, A1, B1]);
MmaTranslator:-Mma:-LeafCount~([A, B, A1, B1]);                   

                    [6209, 120, 3487, 120]

                      [1539, 38, 976, 38]

It will be very difficult (no algorithms available) to simplify A to B.  For the other software too!

eval(A-B, x=2*t):
simplify(factor(expand(%))) assuming real;

                     0

Maple 2024 does not solve the problem.
Here is a workaround:

convert(series((f:=(lhs-rhs)(eq)),t,16),polynom):
solve([coeffs(%,t)]);
eval(f,%); #check

        {A[1]=0, A[2]=0, A[3]=-7/25, A[4]=1/25, A[5]=-2/5, A[6]=1/5, A[7]=3/16, A[8]=0, A[9]=0, A[10]=3/8}
        0

f := n -> (ln(x)^n)^(1/n);

simplify(f(2))  returns  csgn(ln(x))*ln(x)
just because the "simple" function  csgn exists.  csgn(z) equals sqrt(z^2)/z for z<>0.

simplify(f(3))  remains unsimplified because a similar "sign" function does not exist for z^(1/3);

For n = 3 or n=5,  simplify((log(x)^n)^(1/n)) assuming x::positive, x<1   give both  - ln(x) * (-1)^(1/n)
but when n=3,   (-1)^(1/3)   is further simplified to 1/2+I*sqrt(3)/2; 
When n=5, it would be too complicated:  (-1)^(1/5) = sqrt(5)/4 + 1/4 + sqrt(2)*sqrt(5 - sqrt(5))*I/4.

odetest cannot work with such generalized series. But it is easy to test the solution directly: 

fsol:=series(rhs(maple_sol), x=0):
# series(eval(lhs(ode), y(x)=fsol), x);  # O(x^5)
map(series, map(eval, ode, y(x)=fsol), x)

          O(x^5) = 0

Download equal-sloap-curvature-vv.mw