Seems that the code from here will do the trick. The computed values will be stored in the table sols. For the sums from 0 to 5000 you might want to compile R0 and R1. But the solution doesn't seem to change much if I take more terms.

Download fsolve2.mw

Your code computes ContractIndices(g, R0, [[2, 1]]). ContractIndices(R0, g, [[1,1]]) would be

RR[j, k, l, i] := add(R01t[r, j, k, l]*g1t[r, i], r = 1 .. m)

I think something is missing here.

sys := convert(
  [(-x3*cos(x5)+(1010-x4)*sin(x5))^2/x1^2+(-x3*sin(x5)-(1010-x4)*cos(x5))^2/x2^2 = 1,
   ((-50.5-x3)*cos(x5)+(1060.5-x4)*sin(x5))^2/x1^2+((-50.5-x3)*sin(x5)-(1060.5-x4)*cos(x5))^2/x2^2 = 1,
   ((404-x3)*cos(x5)+(1313-x4)*sin(x5))^2/x1^2+((404-x3)*sin(x5)-(1313-x4)*cos(x5))^2/x2^2 = 1,
   -2*x2^2*cos(x5)^2*x3+2020*x2^2*cos(x5)*sin(x5)-2*x2^2*cos(x5)*sin(x5)*x3-2*x1^2*sin(x5)^2*x3-
   2020*x1^2*sin(x5)*cos(x5)+2*x1^2*sin(x5)*cos(x5)*x4 = 0], rational);

Digits := 20:

fsolve(subs(x5 = 2*Pi/5, map(eq -> numer(lhs(eq)-rhs(eq)), sys)), {x1, x2, x3, x4} =~ 0 .. 10^6);
   {x1 = 263.23612356347687652, x2 = 726.46772751454962536, x3 = 583.19908442411991232,
    x4 = 997.50278517973261214}

Perhaps for a fixed x5 there won't be more than one positive solution, but for other values of x5 there will be other solutions.

coulditbe is like an existential quantifier, it answers the question "Do there exist values of the variables that make this predicate true?" If one of the sides of an inequality has a non-zero imaginary part, the inequality will always evaluate to FAIL, not true (even if it's a non-strict inequality). So coulditbe gives false as it should. coulditbe(2+3I<0) also gives false in 2017.2.

These, however, are wrong:

coulditbe(sqrt(x) = -I);
                              true

coulditbe(1/x = 0) assuming x::real;
                              true

There are no values of x for which sqrt(x)=-I and no real x for which 1/x=0 (infinity is realcons but not real).

Regarding the Zeta function, Zeta(s+2I)=0 with real s means there is a root of Zeta(z) with the imaginary part equal to 2. You don't need the Riemann hypothesis to prove that there are no such roots. So s cannot be real.

If it's a system of polynomial equations and inequalities, Maple should be able to solve it:

sys := [y >= 4*x^4+4*x^2*y+1/2, sqrt((1/2)*(x-y)^2-(x-y)^4) = -2*x^2+y^2];

solve([sys[1], sys[2]^2, rhs(sys[2]) >= 0]);
               /            -3\    /           1\ 
              { x = -1, y = -- }, { x = 0, y = - }
               \            2 /    \           2/ 

There are no other real solutions.

Look up Normalizer and Testzero. Without one of them you're not going to get correct results, e.g., Rank(A) will give 54 instead of 53, because Maple won't be able to detect that certain quantities are equal.

In your example it's more efficient to use Normalizer, because evala is able to put the quantities in a kind of a canonical form, so Testzero doesn't need to do anything special:

# with s:=t^3
A, b := GenerateMatrix(sys, var);
A := evala(convert(A, RootOf));

Normalizer := evala;

LinearSolve(A, b); # instantaneous

convert(%, list);
  [_t[1], 0, 0, 0, 0, 0, 3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1]-3*_t[1], 0, 0,
   -3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   -3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1], 0, 0, _t[25], 0, 0,
   3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1]-3*_t[1], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, _t[1], 0, 0, -3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1],
   0, 0, 3*RootOf(_Z^6-_Z^3+1, index = 1)^3*_t[1]-3*_t[1], 0, 0, _t[1]]

Here's how the last equation can be verified. How one would go about deriving it, I have no idea.

Download rho.mw

If you evaluate each line in your SimEquations.mw, it will look like this:

The assignment to eq1 hasn't happened. Evaluate eq1 to see that it's still just 'eq1'.

The reason is that your "eq1:=" is in a separate cell AND it's in a non-executable cell, so you don't get an error when pressing Enter. It's color-coded, but still very easy to confuse it for a cell being executed.

Perhaps you accidentally clicked the Toggle Executable Math button that appears at the top right corner of an input cell.

My suggestion is to use File>New>Worksheet mode, not Document mode, because in Worksheet mode you get the ">" prompts for every input (and  you can still use 2D math).

By the way, the system is homogenous, how can it not have a solution?

I believe you've already asked essentially the same question here, since tanh(x) is trivially related to (exp(4*x)+1)/(exp(4*x)-1).

Maple's symbolic derivative will give the answer for tanh in terms of a sum of Stirling numbers. Pretty much like in your previous question, those sums can be rewritten as a single Zeta(-n) or bernoulli(n+1) or polylog(-n, -1):

(eval(diff(tanh(x), x$n), x = 0))/n!:

subs(Sum((-1)^_k1*factorial(_k1)*Stirling2(n, _k1)/2^_k1, _k1 = 0 .. n) = (2^(n+2)-2)*Zeta(-n), %);
        I*I^(n+1)*2^n*(-1+piecewise(n = 0, 1, 0))*(2^(n+2)-2)*Zeta(-n)*I^n/factorial(n)

The products do not make much sense as written, because you're multiplying monomials of the same form p[i]*x^q[i], so all you can ever get is product(p[i])*x^sum(q[i]).

What exactly are you calculating? g3 has a singularity at 0, so your a is undefined.

If you want to compute the regular part of the Laurent series, I can immediately say that it's g3 minus the singular part, so hello equals g3-1/x :).

If you want the coefficients of the series, it can be done with a bit of cheating:

eval(diff(-2*exp(4*k*x)-2*exp((4*(k+1))*x), x$n), x = 0);
                            n              n
                    -2 (4 k)  - 2 (4 k + 4) 

ga := sum(%, k = 0 .. infinity, formal);
                           (2 n + 2)         
                   ga := -2          Zeta(-n)

so g3=1/x+sum(ga*x^n/n!, n=1..infinity) (notice that here the summation starts at n=1).

In fact, this can be done by hand as well once you look up the generating function for bernoulli(n).

The radius of convergence:

simplify(limit((ga/n!)^(-1/n), n = infinity));
                              1   
                              - Pi
                              2   

What's rather annoying is that evalf cannot compute the sum numerically near x=Pi/2.