My original Answer was wrong because I made M^*.M a symmetric matrix when actually it's hermitian. I thank dharr for pointing out my error. Here's a corrected Answer. The incorrect one is at the bottom of this Answer for the sake of continuity of reading.

I will use L for your l because the default sans-serif font on MaplePrimes essentially makes l unreadable. The trick to doing your problem is expressing L as Lx + I*Ly where Lx and Ly are assumed real, and simplifying conjugated expressions with simplify(..., conjugate) and evalc.

In Maple, M^%H (or M^*) is the HermitianTranspose of M, and A.B is the matrix product of matrices A and B. Spelling out MatrixMatrixMultiply or HermitianTranspose is unnecessary.

Because I use evalc (which assumes variables are real), the assume(a::complex) is needed.

The final eigenvectors below can be further simplified, but it requires subtlety.

Download SymbEigenvectors.mw

Original incorrect worksheet:

Below, I have put your data in exactly the same order as you gave it; however, I've grouped it with square brackets and commas. Please note where the brackets are; the pattern should be obvious. Line breaks, spacing, and indentation are irrelevant; they're only used to enhance readability.

Pts:= [
    [
        [80, 1, 1.80074340195803263], [80, 2, 1.79160585795803287], [80, 3, 1.7805446209580329], 
        [80, 4, 1.76755968995803269], [80, 5, 1.75265106495803280], [80, 6, 1.73581874695803272],
        [80, 7, 1.71706273595803272]
    ],
    [                        
        [85, 1, 1.71097991748770582], [85, 2, 1.70184237348770577], [85, 3, 1.69078113648770581], 
        [85, 4, 1.67779620548770559], [85, 5, 1.66288758048770571], [85, 6, 1.64605526248770562],
        [85, 7, 1.62729925148770591]
    ],                          
    [
       [95, 1, 1.53550477309879818], [95, 2, 1.52636722910871089], [95, 3, 1.51530599205435690], 
       [95, 4, 1.50232106105435669], [95, 5, 1.48741243605435680], [95, 6, 1.47058011805435672], 
       [95, 7, 1.45182410705435672]
    ],
    [                    
       [100, 1, 1.44826461660029512], [100, 2, 1.43912707260029623], [100, 3, 1.42806583560029569],               [100, 4, 1.41508090460029576], [100, 5, 1.40017227960029588], [100, 6, 1.38333996160029579],       
       [100, 7, 1.36458395060029579]
    ]
]: 
plots:-display(PLOT3D(MESH(Pts)), labels= [x,y,z]);

Note that PLOT3D and MESH are all uppercase letters. Some help for them can be found at ?plot,structure; however, the best way to learn about them is to deconstruct the results of plot3d (lowercase) commands with op.

An equivalent command, producing exactly the same plot, is

plots:-surfdata(Pts, labels= [x,y,z])

Either way, the trickier part is gettiing the data arranged correctly in a list of lists of lists (as shown above) or a 3D-Array (as shown by acer). Once that's done, displaying it is easy.

The trick is to insert a fake symbolic exponent and substitute 1 for it after the extraction, like this:

Powers:= (e, T::type)-> local One;
    eval(
        indets(combine(evalindets(e, T, `^`, One), power), T^anything),
        One= 1
    )
:
e:=_C1^2+_C1+_C1*_C2^3+a+b:
Powers(e, suffixed(_C, nonnegint));
                       /        2     3\ 
                      { _C1, _C1 , _C2  }
                       \               / 

By the way, suffixed(_C, ...) is a subtype of symbol, so your And(symbol, ...) is redundant.

It's not a caching issue. What's happening is that the elapsed time between your calls to randomize is much shorter than the "ticks" of the clock. There's never a need to call randomize more than once. Indeed, the effect of calling randomize twice in a short time period is to reset the random-number generator to a specific entry in its sequence, so that's about as unrandomized as you can get.

My feeling about randomize is that it should only be called at the top level, never in a procedure.

Also, there's never a need to pass randomize a custom clock value. The simple call randomize() (with no arguments) already uses a clock-based value that's sufficient for any purpose.

Like this:

collect(T, indets(T, function));
                    f(x) + 3 + (b + 1) g(x)

@Kitonum is correct that shading regions of polar plots is difficult in Maple, and a "polygon" approach is needed. Here's a way to do it that's perhaps more intuitive than his: I extract the polygon vertices directly from the output of polarplot commands. There's no need to look at the numerical vertices; it's all done rather automatically. Some additional benefits of my approach are that I put the plots on Maple's polar axes rather than Cartesian axes, and no explicit coordinate conversions are required.

I also show three distinct ways to get the area.

Download PolarArea.mw

Use sort with this key:

sort(F, key= [f-> subs(0= infinity, degree~(f, [x4,x3,x2,x1])), nops])

I assume that by "degree of nonlinearity in x4" you mean nothing more than simply "degree in x4 (regardless of the presence of other variables)". If that's not what you mean, let me know. Also, this solution assumes that the polynomials are expanded.

Change
yB_:= y -> sol4[1]:
to
yB_:= MakeFunction(sol4[1], y):

You are correct that the y in the solution is not being seen. This is because there is no y explicitly on the right side of the arrow. In order to make a non-explicit variable into a procedure parameter, you need to use MakeFunction instead of ->.

Creating the set of equations can be done as simply as this:

{seq}(Matrix(2, symbol= a) - Matrix(2, symbol= b) =~ Matrix(2, shape= identity));
        {a[1, 1] - b[1, 1] = 1, a[1, 2] - b[1, 2] = 0, a[2, 1] - b[2, 1] = 0, a[2, 2] - b[2, 2] = 1}

@michaelvio 

The vast majority of the slowness that you were experiencing was due to your use of int (symbolic integration) versus Int (numeric integration, when it's used in the context that you were using).

It's true that using a larger epsilon and a smaller numpoints will save some time, but that savings is infinitesimal in this case. I recommend that you remove both options so that those things will be done at their default values (although epsilon= 1e-6 is sufficient for plotting purposes).

You should also set Digits:= 15 right after the restart at the top of the worksheet. The extra precision makes a huge difference in the final plot, and the extra time needed isn't even noticeable.

To construct the function J1 as you just described, all that you need to do is

J1:= unapply(subs(x0= 1/t0, J), t0);

This is your last plot at the default Digits:= 10:

And this is using Digits:= 15:

It's a similar issue as with your recent Question regarding applying eval to Int or Intat: eval uses rules, usually to ensure the mathematical validity of the proposed substitution. The rule in this case is the procedure `eval/limit`. To create a more-sophisticated rule, read the procedure (showstat(`eval/limit`), it's a fairly easy read with only a moderate amount of Maple knowledge, and it's only 42 lines (Maple 2023)) and either write an overload or, perhaps, modify the procedure directly (ToInert, FromInert).

If you'd like it to follow rules additional to what's already in the procedure, and you can state them for me with a moderate amount of mathematical precision and generality, I may be able to suggest a way to code it. 

Let's consider your examples individually:

1. e:= A+B*limit(f(x), x= infinity);  eval(e, limit(f(x), x= infinity)= 1);
This case doesn't involve changing anything inside any limit. Thus, the rule isn't used (verified by using trace(`eval/limit`)), and you got your expected result.

2. e:= A+B*limit(2*f(x), x= infinity);  eval(e, limit(f(x), x= infinity)= 1);
In this case you're asking it to evaluate an expression containing a limit by using information about another limit (albeit, a very closely related one); the eval rule will be used. I think that the mathematical rule (i.e., abstract, not Maple) that you're implicitly invoking here is always valid, although I'm not 100% sure about all weird "corner cases", domains other than 2-sided real, and multivariate (more than one bound variable (e.g., limit(..., (x,y)= ...)) cases. That rule, stated generally, is that factors that don't depend on the bound variable(s) can be factored out of the limit. Command expand will do this, so a quick solution is 
eval(expand(e), limit(f(x), x= infinity)= 1)

[to be continued]

The correct command is is: 

is(sqrt(3) < 3);

Presumably it just expands things to a sufficient number of decimal places to make sure.

You wrote:

• Should one then use is(n>=m) instead of evalb(n>=m) always?

No! That's inefficient. evalb is builtin (kernel code rather than Maple). is is a highly symbolic command written in Maple. In the vast majority of cases, you know a priori that both sides are type numeric, and is isn't needed.

• I can also always apply evalf() on each side before. But this seems like an overkill to me.

I think that's the approach that is takes, and it might increase Digits if the user's set value is not high enough to make absolutely sure of the inequality.

I said that I could probably come up with something better for large n and m < < n, and here it is.

To make the algorithm fast, there are three things to avoid:

1. finding all the large prime factors of n when that is difficult to do,
2. finding all the divisors and then filtering them when the number of divisors is large,
3. checking every positive integer < m when m is large.

The command ifactors(n, easy) finds only the prime factors that are easy to find.

Download SmallDivisors.mw

This may be suitable for you:

select(`<`, Divisors(n), m)

Surely it's more efficient than what you were doing. For an easy-to-factor n >> 2^64 and m << sqrt(n), I could probably come up with something even better.

My guess is that you actually want the simplified product of an arbitrary number of factors of T[r]*P[r] --

T[1]*P[1]*...*T[n]*P[n]

-- because otherwise your Question would be exceedingly trivial. But first, you should stop trying to use an index (your r) as a parameter. Maple does allow the use of indices as parameters, but it's a somewhat more advanced topic than your current level of Maple knowledge. So, for now, let's just make r an ordinary parameter:

T:= r-> exp(r+1);  P:= r-> exp(r-1);

Your product can be represented as an inert product with the Product command:

Product(T(r)*P(r), r= 1..3);

Now, here's the non-trivial (but still conceptually easy) part: Make the upper bound, 3, arbitrary, and make the product non-inert by using product (lowercase):

product(T(r)*P(r), r= 1..n);

Of course that can be greatly simplified:

simplify(%);

By the way, can anybody give me an English word that's actually used in mathematical prose that completes this formal two-level analogy?:

(Addition is to summation as multiplication is to __?__)
as in Maple
(add is to sum as mul is to product).

That word, if it exists, would replace "products" in this Answer's Title.