@Carl Love 

A "Godelian" example would look quite different. This one should be easy for a human but difficult for Maple. BTW, the answer is "true".

A tough test for Maple would be:

F := piecewise(exp(x+exp(x)-sin(x^2))>=x+1-cosh(x/2), 2*x+1, x^2);

@Kitonum 

You are right. I was fooled by Maple's solve. I am disappointed that Maple fails for a simple rational polynomial system:

sys := [x^2 + y^2 - x*y - 1 = 0, y^2 + z^2 - y*z - a^2 = 0, z^2 + x^2 - x*z - b^2 = 0]:
ab:={a=10,b=104/10}:
solve(eval(sys,ab), allsolutions):
evalf[20]([%]);

    [{x = 1.1419533726022634876, y = 0.42279817537010308072, z = -9.7818952278636396166}]

The option size for 3D plots (implemented via ROOT) does not work as expected (unlike the 2D case). Here the plot is cropped sometimes; the PLOT3D driver is probably not yet ready for this option.
But of course it's nice to have it now, even if we must play with the size parameters in order to obtain an acceptable plot.
Examples:

Compare:
plot3d(x^2+y^2, x=-2..3,y=-2..3);
plot3d(x^2+y^2, x=-2..3,y=-2..3, scaling=constrained);

with

plot3d(x^2+y^2, x=-2..3,y=-2..3, size=[400,800]);

plot3d(x^2+y^2, x=-2..3,y=-2..3, scaling=constrained, size=[400,800]);

plot3d(x^2+y^2, x=-2..3,y=-2..3, scaling=constrained, size=[300,1800]);

@amrramadaneg 

I don't know what you have computed in Mathematica, but:

1. U1(y,z)  does not depend on y.
2. U1 exists only vor small z. For example, for z=0.5, the integral does not exist.
plot(U(x,0.5), x=0..1, view=-5..5);

This is a mathematical fact, not related to any CAS.

@amrramadaneg 

Your plot does NOT make sense because int without an interval as parameter will introduce arbitrary constants.
Kitonum changed this by taking definite integrals using the intervals 0.1 .. X,  0.1 .. Y.
Actually 0.1  is not necessary because x=0 is a removable singularity.
The problem is the next singularity at x = x0 = 0.42...  For x>=x0  the integral does not exist, so your x = 0..10 (in plot) is nonsense too. The computations takes long because Maple tries to compute a lot of divergent integrals.

@Kitonum It does not answer the question.
For m=0 the problem is trivial. We have m=1.

restart;
S:=msolve(x^3+y^2 = 123, 124):
S[rand(1..nops([S]))()];
                        {x = 83, y = 72}
eval(x+y, %);
                              155

@Stretto 

Now, that you know the syntax, you have options.
For example, in Maple 2019:

(for i from 1 to 3 do    
    A := i:
    B := 10*i;
end);

        10, 20, 30

Or, you may redirect the output into a file and then filter its content.

@Carl. I was intrigued by the results obtained by your procedure and I have looked closer to the u[i] inequalities.
The conclusion is that they are ok not only for n, but also for any n'>=n. So, your LPSolve call is correct and should produce the solution. Congratulations for your intuition and vote-up!
Unfortunately, LPSolve does not obtain the global minimum. I have checked with my brute-force procedure; it gives a better solution but this solution satisfies the LPSolve conditions, so it should have been found!

@Carl Love 

1. The triangle inequality must be checked for all i,j,k,  not only for i<j<k. In your example, the resulting matrix violates the triangle inequality.

2. The u[i] variables are for the TSP problem, not VRP. There, u[i] need not be >=0 [actually, obviously, if u[i] satisfy the conditions, u[i]+const satisfy too]. A special solution for u[i] is the step at which the vertex i is visited (but u[i] is not unique).
In VRP the problem is that the tours do not have the length n, so, insted of n should be the lenth (which is not known, probably it is possible to use the sum over i,j of x[i,j,v]).
There is another problem. In my brute-force solution, the tours are disjoint (except vertex 1, of course). I don't have an example where nondisjoint tours lead to a better solution.

@Carl Love 

When running your worksheet I obtain other results. Please check.
Note also that " Check and ENFORCE triangle inequality " does not work this way.

@Carl Love You forgot op:
factor(op(1, convert(cos(Pi/257), RootOf)), sqrt(257));

@Kitonum 

Sure.
Note that in fact (-1)^(1/7) is itself a radical, so, if we want a simple cos(Pi/7)  we may use ( (-1)^(1/7) + 1/(-1)^(1/7) )/2.

@Kitonum 

Actually,