The Maple implementation NumberTheory:-Totient is of course faster (for large numbers), but the simplest one is probably:

eulerphi:=(n::posint) -> add(`if`(igcd(k,n)=1,1,0),k=1..n);

Should be

solve(SysEqE,{diff(phi[1](t),t),diff(phi[3](t),t)});

op(0,a);

But you should try to see how this appeared, because it is a nonsense (probably from a programming error).

I don't think it has a special name; it is a multiplication by a diagonal matrix.

But are you sure it is useful? AFAIK the eigenvector algorithms always do a row (or column) normalization.

Use
evalc(Re(expr));

Download sysode-int.mw

This is a maths problem, not a Maple one. See https://en.wikipedia.org/wiki/Cycles_and_fixed_points

Of course it is easy to implement any of the formulae.

Here is a direct solution (based on simplex) asked by the OP who has an older version of Maple.

FacetsOfPolyhedralCone:=proc(L::list)
local i,j,J,A,X,u;
X:=indets(L,name);
A:=LinearAlgebra:-GenerateMatrix(L,X)[1]:
J:={seq(1..nops(L))}:
for j in J do
  add(A[j]-u[i]*A[i],i=J minus {j}); convert(%,list)=~0;
  if simplex[minimize](0, %, NONNEGATIVE) <>{} then J:=J minus {j} fi;
od:
J, L[[J[]]];
end:

L := [y-z, 3*y-2*z, 2*y-2*z, x-2*y+z, x-y, 2*x-y, x-z, x+y-z, x, y, z]:  # >= 0

FacetsOfPolyhedralCone(L);
       {3, 4, 11}, [2*y-2*z, x-2*y+z, z]

solve({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])}, {x[1], x[2], x[3], x[4]}); 

         {x[1] = 0, x[2] = 0, x[3] = 0, x[4] = 0}

It is a bug. The result of eliminate should be:

Mathematically it should be also added: x[5]<>0, x[7]<>0, x[4]<>0,  x[6] <> -(2*x[7]+1)/(x[5]*x[7])

@mmcdara 

Let's do the computations for 3 digits.

restart;
Digits:=3:
a,b,c:= sqrt~([2,3,6])[]:
fa,fb,fc := evalf(%);
                 fa, fb, fc := 1.41, 1.73, 2.45
fa*fb;
                              2.44
# Do you expect  2.45?  How?
141*173;
                             24393

num1dsubspaces := (p,n) -> (p^n - 1)/(p-1)

ineqs:=
y-z >=0,
3*y-2*z >=0,
2*y-2*z >=0,
x-2*y+z >=0,
x-y >=0,
2*x-y >=0,
x-z >=0,
x+y-z >=0,
x>=0,
y>=0,
z>=0:
with(PolyhedralSets):
A := PolyhedralSet([ineqs]);
B := PolyhedralSet([ineqs[3..4],ineqs[-3..-1]]);

evalb(convert(A,string)=convert(B,string));

        true

Add:

f1,f2,f3:=seq(unapply('evalf'(x__||i),t2), i=1..3):
plot([f1,f2,f3], 0..0.2);

You will not be able to plot over  0 .. Pi/2  because the integrals are not convergent here.