It does seem to be a bug. When using RootOf without a selector, there is some ambiguity as to how that object is interpreted (the help vaguely says that "Some functions use only the principal root"). But here each RootOf is uniquely defined, and still the same problem occurs:

LinearAlgebra:-Eigenvectors(<0, RootOf(_Z^2-2, index = 1); RootOf(_Z^2-2, index = 2), 0>);
        Vector[column](2, [RootOf(_Z^2+2, index = 1), RootOf(_Z^2+2, index = 2)]),
        Matrix(2, 2, [[0, 0], [0, 0]])

If you want an answer in terms of your auxiliary polynomials, you can use Groebner:-NormalForm:

rels := [Ones-(i2+i3+i4+i5), Twos-twos, Threes-threes, Infosquare-infosquare, Erbij-erbij]:
bas := tdeg(i2, i3, i4, i5, Ones, Twos, Threes, Infosquare, Erbij, i1):
gb := Groebner:-Basis(rels, bas):

Groebner:-NormalForm(noemer, gb, bas);
        -Ones*Threes*i1+Twos*i1^2+Erbij*Ones-2*Erbij*i1+Threes*Infosquare

Groebner:-NormalForm(Teller, gb, bas);
        -2*Ones*Threes*i1^2+2*Twos*i1^3+2*Erbij*Ones*i1-4*Erbij*i1^2-4*Ones*Threes*i1+
        Ones*i1*Infosquare+2*Threes*i1*Infosquare+4*Twos*i1^2+i1^2*Infosquare+4*Erbij*Ones-
        8*Erbij*i1+(4/3)*Ones*Threes+4*Threes*Infosquare-(4/3)*Twos*i1+4*Erbij

This is a bit of a mess. But first, some preliminaries. Consider:

sum(rand(0 .. 1)(), k = 1 .. 20);
                               20

add(rand(0 .. 1)(), k = 1 .. 20);
                               9

That is, even if the iterator 1..20 is numeric, sum evaluates the expression once for symbolic k; add doesn't. Also,

0^k;
                               0
0^0;
                               1

So if your summand contains 0^k, evaluating it for symbolic k will be different from evaluating it only for numeric values of k. And when you use unapply, you evaluate the sum from 1 to 5 for symbolic x, x^0 becomes 1, and that 1 won't disappear regardless of what you substitute for x later. So far, so good. But:

0^(k-1);
                            0^(k-1)

simplify(%);
Error, (in simplify/commonpow) numeric exception: division by zero

which isn't consistent with 0^k. So I'm actually not sure how sum simplified your f(x) to 0. Further,

restart;
 
(x -> sum(x^(k-1), k = 1 .. 5))(0); # in Maple 2017.3
                             0

sum(0^(k-1), k = 1 .. 5);
                             1

(x -> sum(x^(k-1), k = 1 .. 5))(0); # identical to the first input
                             1

The second input makes the two other inputs, which are identical, give different results.

There is ToInert, but the output will be pretty verbose. Here is an alternative:

inert := proc(e) local head, params;
  if e::atomic then return e end;
  params := op(map(inert, e));
  head := inert(op(0, e));
  (h -> 'h'(params))(eval(head, 1)) # put head into uneval quotes
end:

inert('seq([i], i = 1 .. 2)');
                     seq(list(i), `=`(i, ..(1, 2)))

inert('(f+rand())(x/rand())');
               `+`(f, rand())(`*`(x, `^`(rand(), -1)))

inert('(e->e)(x)'); # eval(head, 1) in the code is needed for this case
                          (e -> e)(x)

There are some issues (like the range operator being shown as .. instead of `..`, which will not be interpreted correctly on input), but mostly that seems to do the job.

Then you can write some code to draw a graph with those nodes.

Seems to be just a question of numeric precision. You can do your computation with exact quantities (Ec:=convert(0.380e12, rational, exact), etc.), then for m=10, the determinant of MT will be a polynomial of degree 8 with rational coefficients, and you don't need to increase Digits to compute its roots. While if you start with appoximate values, you need a high level of precision to get a sufficiently close polynomial.

Another thing is the precision of the input values. If your Ec can be 0.3800e12 or 0.3804e12, and the corresponding difference in the final result is 0.005, then it doesn't make much sense to compute the next digits.

You need to compute the bordered Hessian, which is just the Hessian of the Lagrange function wrt all its variables:

Student:-VectorCalculus:-Hessian(f(x,y,z)+l[1]*(x*y-1)+l[2]*(y^2+z^2-1), [l[1], l[2], x, y, z]);

but the conditions for determining the type of an extremum are slightly different from the unconstrained case. Here the value of the determinant will be either +4 or -4 at the extrema. The test says that the points where it's -4 are local maxima, and the points where it's +4 are local minima.

There are two local maxima (since the problem is symmetric about the origin) and two local minima.

For numeric values of alpha and mu, you can compute the Fourier coefficients directly:

int(eval(Ia*exp(-I*k*theta), [alpha = 524/1000, mu = 549/1000]), theta = 0 .. 2*Pi);

The output is fairly large because Ia[Int1],...,Ia[Int13] are undefined in your worksheet, but it appears to be correct.

If you write your roots of unity as exp(I*a) instead of just a, it'll make the simplifications easier:

ubound := proc(polyq, q) local cs, ms, i;
  cs := Vector([coeffs](polyq, q, ms));
  for i to LinearAlgebra:-Dimension(cs) do
    cs[i] := expand(cs[i]);
    if cs[i]::`+` then cs[i] := [op](cs[i]) else cs[i] := [cs[i]] end;
    cs[i] := LinearAlgebra:-Norm(Vector(cs[i]), 1)
  end;
  cs . Vector([ms])
end:

ubound(q^3-exp(I*a)*q^2+exp(I*b)*q^2-5, q);
              5+q^3+(exp(-Im(a))+exp(-Im(b)))*q^2

Then evalc(%) (which by default assumes that the variables are real) or simplify(%) assuming a::real, b::real will do the job.

expr732vc and expr732wc are equations (op(0, expr732vc) is `=`), coeff can't be applied to equations.

(lhs-rhs)(expr732vc) gives 0., the same with expr732wc. I don't know if that's intended.

Maybe you want coeff~(lhs(expr732vc)-lhs(expr732wc), epsilon, [0, 1]).

The easiest way is probably to work with vectors and use the parametric plot:

scale := v -> convert(Matrix([[2, 1], [1, 2]]).Vector(v), list):

proc() local rng := -5 .. 5;
  plots:-display(
    plot([seq](op([[t, i, t = rng], [i, t, t = rng]]), i = rng), linestyle = dash, color = gray),
    plot([seq](op([[(op@scale)([t, i]), t = rng], [(op@scale)([i, t]), t = rng]]), i = rng), color = gray),
    plot([[t, t/2, t = rng], [(op@scale)([t, t/2]), t = rng]], thickness = 2, legend = ["original", "scaled"]),
    view = [rng, rng])
end();

The equation for x looks like this:

elim := eliminate(sys, {H, y});
      [{H = 131+104*ln(x-40), y = -126000000/4109923+x},
      {786000+624000*ln(x-40)+4181*ln((1/4109923)*(-290396920+4109923*x)/(x-40))}]

This is a sum of two logarithms plus a positive constant, thus the argument of the second logarithm needs to be small to compensate for the positive terms. Search for the solution in the form x0+eps:

series(eval(elim[2, 1], x = 290396920/4109923+eps), eps = 0, 1);

solve(convert(%, polynom), ln(eps));
      -786000/4181-(624000/4181)*ln(126000000/4109923)-ln(4109923/126000000)

The relative error in eps will also be on the scale of eps, so the values of x,y,H will have about 600 correct digits.

The general way is to compute the Groebner basis with the specified elimination order.

Ignoring the denominators, we obtain:

Groebner:-Basis((numer@(lhs-rhs))~(
  {x[5] = x[2]/x[1], x[6] = x[3]/x[2], x[7] = x[1]/x[4], x[8] = (2*x[2]+x[4])/(2*x[1]+x[3]+x[4])}),
  lexdeg(index~(x, [$1..4, 8]), index~(x, [$5..7])));
                  [-x[4]*x[5]*x[6]*x[7]+x[3], -x[4]*x[5]*x[7]+x[2], -x[7]*x[4]+x[1],
                  x[4]*x[5]*x[6]*x[7]*x[8]-2*x[4]*x[5]*x[7]+2*x[4]*x[7]*x[8]+x[4]*x[8]-x[4]]

eliminate(), on the other hand, doesn't seem to be a well-defined operation (e.g., I still don't know whether the second eliminate example here works as expected or not).

Take a simpler example:

series(sqrt(a+b*u^2), u = 0, 3);
                 sqrt(a)+(1/2)*b*u^2/sqrt(a)+O(u^4)

There is no linear term -- as long as a is not zero. In nu3, instead of a, there is an expression that is zero only for some values of w. You can add the corresponding assumption (which is valid for w>sqrt(2)-1):

series(nu3, u = 0, 3) assuming -2*sqrt(w)*sqrt(w+2)+2*sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)+2 = 0;
      -sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)+w-1+sqrt(2)*sqrt((2*sqrt(w)*sqrt(w+2)*sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+
      2*w+1)+2*sqrt(w)*sqrt(w+2)-2*w*sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)-sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)-
      2*w-1)/sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1))*csgn(sqrt((2*sqrt(w)*sqrt(w+2)*sqrt(-2*sqrt(w)*sqrt(w+2)+
      w^2+2*w+1)+2*sqrt(w)*sqrt(w+2)-2*w*sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)-sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+
      2*w+1)-2*w-1)/sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1))*u)*u-(1/2)*(4*sqrt(w)*sqrt(w+2)-4*w-2)*u^2/
      sqrt(-2*sqrt(w)*sqrt(w+2)+w^2+2*w+1)+O(u^3)

Now it has the correct linear term (times csgn(u)).

A cleaner way is to rewrite nu3 as a RootOf:

nu3r := RootOf(evala(Norm(_Z-nu3)));

series(nu3r, u = 0, 3);
      w-sqrt(w^2+2*w)+RootOf(_Z^2*sqrt(w^2+2*w)-_Z^2+4*w*sqrt(w^2+2*w)-4*w^2+2*sqrt(w^2+2*w)-8*w)*u-
     (2*sqrt(w^2+2*w)-2*w-1)*u^2/(sqrt(w^2+2*w)-1)+O(u^3)

and you just need to select the correct RootOf branch.

Adding expand around binomial is faster, everything is polynomial in a and b right away:

det := (n) -> LinearAlgebra:-Determinant(Matrix(n, n, (i, j) -> expand(binomial(a*i+b*j, j))));

det2 := (n) -> a^(n*(n-1)/2)/n!*(n!*a^n+b*add(abs(Stirling1(n+1, i))*a^(n+1-i), i = 2 .. n+1));

expand([seq](det(i)-det2(i), i = 1 .. 10));
                 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

The closed form is

det3 := (n) -> a^(n*(n+1)/2+1)/n!*(n!/a+b*((1/a+n)!/(1/a-1)!-n!/a));

(expand@@2)([seq](det(i)-det3(i), i = 1 .. 10));
                 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Groebner:-NormalForm gives the form you want:

eu := -2*n^2*B^2+2*n*B^2*p-2*n*B^2*p^2+n*B^2-B^2*p+B^2*p^3+n*c*d;

Groebner:-NormalForm(eu, [n-p+p^3, n-p+p^2], tdeg(n, p), 'quos');
                             n*c*d
quos;
                        [B^2, -2*n*B^2]

But the divisors do not form a basis, so the remainder is not unique (same as in this question).

Hope that it will not only simplify but also expand your life.