The first integral is calculated incorrectly (a bug). If test1 to split into two integrals, then all is well:

evalf(Int(cot(2*t-1-I), t=0..0.5))+evalf(Int(cot(2*t-1-I), t=0.5..1));

                                         0.+1.115976093*I

You can generate a random Hermitian matrix as follows:

LinearAlgebra[RandomMatrix](4, generator=rand(-9..9)+I*rand(-9..9), outputoptions=[shape=hermitian]);

A:=x/2+y/2:

subs(a=``, normal(A/a));

I know only the way by  textplot  command:

test := proc(x)

local A, B;

A:=plots[textplot]([0,0, `x is equal to:  `], color=green, font=[TIMES, ROMAN, 18], align=left);

B:=plots[textplot]([0,0, 2], color=black, font=[TIMES, ROMAN, 18], align=right);

plots[display](A, B, axes=none);

end proc:

Example:

test(2);

test := proc (x)

cat(`x is equal to:  `, x);

end proc:

Example:

test(2);

                x is equal to:  2

A:=proc(M)

Matrix(M, M, {(1,1)=1/2, (1,2)=1/2/sqrt(2), (2,1)=-1/8/sqrt(2), seq((i,1)=-1/2/sqrt(2)/i/(i-2), i=3..M), seq((i,i-1)=-1/4/(i-2), i=3..M), seq((i,i+1)=1/4/i, i=2..M-1)});

end proc:

Examples:

A(3), A(4), A(5);

If the parameter values are specified​​, the equation maybe easily solved numerically by finding all the roots. To simplify the typing I changed the notations of variables and parameters. At the output all are unchanged.

Example:

x:=alpha: t0:=mu0: t1:=mu1: s:=sigma:

eq:=60*t0*x^2*s^2-60*x*s^2*t0^2+60*t1*s^2*x^2-60*s^2*t1^2*x-15*t0*x^2-15*t1^2*x^4+30*t1^2*x^3-30*t1^3*x^2+15*t1^4*x-40*x^3*s^2+30*x^3*t0^2-30*t0^3*x^2+240*x*s^4+15*x*t0^4+6*x^5+20*s^2-60*t1*s^2+60*t1^2*s^2-120*s^4+15*t1-30*t1^2+30*t1^3-15*t1^4-3=0; 

fsolve(eval(eq, {mu0=1, mu1=2, sigma=3}), alpha, complex);

The easiest way - translation of all constants to symbols and use the concatenation operator.

Examples:

`3`+`7`+`1`=`11`;

`3+7+1`=`11`;

`(-2)`||` (-1)`/`(-5)`||` (-1)`;

add(a[i]*mul(s-b[j], j={$1..n} minus {i}), i= 1..n);

Example:

n:=5:
add(a[i]*mul(s-b[j], j={$1..n} minus {i}), i= 1..n);

PS:  In Standard instead of  minus  you can use  \  :

n := 5:  add(a[i]*mul(s-b[j], j={$1..n}\{i}), i=1..n);

There are the functions  max  and  min  but not  Max  and  Min.  

The equation  max(min(x, 1), min(x, 2))=a  is equivalent to  

min(x,1)<=a  and  min(x,2)<=a  and  (min(x,1)=a  or  min(x,2)=a)

so

solve(min(x,1)<=a  and  min(x,2)<=a  and  (min(x,1)=a  or  min(x,2)=a), x);

eq is your last equation in which all terms are moved to one side.

eq := 2*sin(alpha(t))*cos(a[1])*h[1]*rBTP[1]-2*cos(alpha(t))*cos(beta(t))*h*h[1]+2*sin(beta(t))*sin(a[1])*h*rF[1]+2*cos(beta(t))*cos(gamma[1](t))*e[1]*rBTP[1]-2*cos(alpha(t))*cos(a[1])^2*rBTP[1]*rF[1]-2*z(t)*sin(alpha(t))*cos(a[1])*rBTP[1]+2*z(t)*cos(alpha(t))*cos(beta(t))*h+2*cos(beta(t))*cos(a[1])^2*rBTP[1]*rF[1]+2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*cos(gamma[1](t))*e[1]*rBTP[1]-2*z(t)*h[1]+2*sin(gamma[1](t))*e[1]*h[1]-2*z(t)*sin(gamma[1](t))*e[1]-2*cos(beta(t))*rBTP[1]*rF[1]-2*cos(gamma[1](t))*e[1]*rF[1]-l[1]^2+2*sin(alpha(t))*cos(beta(t))*cos(a[1])*cos(gamma[1](t))*h*e[1]-2*cos(alpha(t))*sin(beta(t))*sin(a[1])*sin(gamma[1](t))*e[1]*rBTP[1]-2*sin(alpha(t))*sin(beta(t))*cos(a[1])*sin(a[1])*rBTP[1]*rF[1]-2*sin(beta(t))*sin(a[1])*cos(gamma[1](t))*h*e[1]+2*cos(alpha(t))*cos(a[1])^2*cos(gamma[1](t))*e[1]*rBTP[1]+z(t)^2+h[1]^2+rBTP[1]^2+e[1]^2+rF[1]^2+h^2-2*cos(alpha(t))*sin(beta(t))*sin(a[1])*h[1]*rBTP[1]-2*cos(alpha(t))*cos(beta(t))*sin(gamma[1](t))*h*e[1]-2*cos(beta(t))*cos(a[1])^2*cos(gamma[1](t))*e[1]*rBTP[1]+2*z(t)*cos(alpha(t))*sin(beta(t))*sin(a[1])*rBTP[1]-2*sin(alpha(t))*cos(beta(t))*cos(a[1])*h*rF[1]+2*sin(alpha(t))*cos(a[1])*sin(gamma[1](t))*e[1]*rBTP[1]:
eq1 := collect(subs(sin(gamma[1](t)) = x, cos(gamma[1](t)) = y, eq), [x,y]);
A := coeff(eq1, x):
B := coeff(eq1, y):
C := tcoeff(eq1, [x, y]):
sol := solve(a*sin(t)+b*cos(t)+c, t):
gamma1 := unapply(eval(sol[1], {a = A, b = B, c = C}), t);  # 2 series of solutions
gamma2 := unapply(eval(sol[2], {a = A, b = B, c = C}), t);

It would be better to use  solve  instead of  isolate, because  isolate does not provide all solutions . See the example:

eq:=a*sin(t)+b*cos(t)+c=0:

isolate(eq, t);

solve(eq, t, AllSolutions);

This is a cylindrical body, in the base of which the circle  x^2+y^2=2*x  or  (x-1)^2+y^2=1, and on the top  it is limited by a paraboloid of revolution x^2+y^2=z .

plot3d(x^2+y^2, y = -sqrt(2*x-x^2) .. sqrt(2*x-x^2), x = 0 .. 2, scaling = constrained, filled = true, axes = normal, numpoints = 10000);
Volume=int(x^2+y^2, [y = -sqrt(2*x-x^2) .. sqrt(2*x-x^2), x = 0 .. 2]);

restart;

X:=[3,4,5,6,7]:

Y:=[7.42494922444550, 3.67768248674133, 2.52235142453921, 1.95610223891559, 1.61770309810016]:

A:=plot(X,Y, style=point, symbolsize=20):

CurveFitting[LeastSquares](X,Y, x, curve=a/(x-2)^0.94);

B:=plot(%, x=0..8, -1..9, discont=true,color=blue):

plots[display](A,B, scaling=constrained);  

PS.  Your attached file is blank. 

I was wrong about "a lot of time." Still, the capabilities of modern computers boggle my mind ! Less than 7 minutes the program went through  21^6 = 85766121 variants and  found all solutions. The first 10 and last 10 solutions displayed.

restart;

ts := time(): N := 0:

for c from -10 to 10 do
for d from -10 to 10 do
for e from -10 to 10 do
for f from -10 to 10 do
for g from -10 to 10 do
for h from -10 to 10 do
a := c+2*d+4*e+6*f+10*g+22*h-16: b := -2*c-3*d-5*e-7*f-11*g-23*h+20:
if a >= -10 and a <= 10 and b >= -10 and b <= 10 then N := N+1:
Sol[N] := [a, b, c, d, e, f, g, h] end if:
end do: end do: end do: end do: end do: end do:

N; seq(Sol[i], i = 1 .. 10); seq(Sol[i], i = N-9 .. N);
time()-ts;