@Bendesarts  The equation  eq_liaison  can be simplified by one-lined code:

eq_liaison:=sin(gamma0(t))*cos(beta0(t)) = -(sin(gamma[1](t))*sin(psi[1](t))*sin(theta[1](t))-sin(gamma[1](t))*cos(theta[1](t))*cos(psi[1](t))+cos(gamma[1](t))*sin(psi[1](t))*cos(theta[1](t))+cos(gamma[1](t))*cos(psi[1](t))*sin(theta[1](t)))*cos(beta[1](t));

lhs(eq_liaison)=1/2*(trigsubs(2*combine(rhs(eq_liaison))))[];

I do not understand your desire maximally to automate the task of fairly particular. I believe that the most effective way to work in Maple - is a reasonable combination of manual and automatic operation. The criteria of effectiveness - is the most quick and easy solution to a specific problem.

@Carl Love  Good solution! Unfortunately, I do not know much English and for me often easier to write new code than to understand the capabilities of the existing built-in commands.

@vv  I agree, but I understand the OP's question as finding the exact value  N(epsilon)

@Age  The quotes in  'deriv'(x)  exclude the premature  calculation. Without quotes  x  in  deriv(x)  should be a constant.

See

deriv(x);

Error, (in deriv) invalid input: fdiff expects its 2nd argument, N, to be of type {integer, name, list(integer), list(name), list(name = constant), set(name), set(name = constant), name = constant}, but received t = x

@Bendesarts 

plots[sparsematrixplot](Phi, matrixview, tickmarks=[[seq(i=q||i, i=1..3)], [seq(-i=Eq||i, i=1..3)]], scaling=constrained);

@Earl  I do not have Maple 15, but in the standard interface Maple 16 this is working properly.

@Bendesarts   I think that there is no single simplifying commands for all occasions. Your last expression can easily be simplified in a similar manner, applied to the right side:

-cos(gamma0(t))*cos(beta0(t)) = -(sin(psi[1](t))*cos(theta[1](t))*sin(gamma[1](t))-sin(psi[1](t))*sin(theta[1](t))*cos(gamma[1](t))+cos(theta[1](t))*cos(psi[1](t))*cos(gamma[1](t))+cos(psi[1](t))*sin(theta[1](t))*sin(gamma[1](t)))*cos(beta[1](t));

lhs(%)=(trigsubs(-2*combine(rhs(%)))/2)[];

@Leinad  What you asked about:

k:=4+f*I:

evalc(abs(k)-2=0) assuming f::real;

                                            (f^2+16)^(1/2)-2 = 0

@vv  Thank you, I didn't know before about  trigsubs  command.

@Bendesarts  Automatically you can do so:

S:=indets(eq_liaison);

S1:=indets(select(t->has(t, sin) or has(t, cos), S));

S minus S1 ;

@Bendesarts  Substitution  T  replaces more complex names with  simple names. I do not know why without these changes  the simplification is not fully implemented.

@acer  Thank you for your attention! Of course should be  s^(11/2)   instead.

@nMaple  Here is an animation with front moving part filled with yellow color (the initial yellow region bounded by y=x^2+1, x=0, and y=5 and rotation about x=-1):

restart;

with(plots): with(plottools):

ax := spacecurve([-1, t, 0], t = -1 .. 6, color = red, thickness = 2, linestyle = dash):

s := plot3d([x, y, 0], x = 0 .. 2, y = x^2+1 .. 5, style = surface, color = yellow):

F := q->rotate(q, alpha, [[-1, 0, 0], [-1, 1, 0]]):

s0 := animate(display, [F(s)], alpha = 0 .. 2*Pi, frames = 60):

s1 := animate(plot3d, [[(x+1)*cos(t)-1, x^2+1, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

s2 := animate(plot3d, [[(x+1)*cos(t)-1, 5, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

s3 := animate(plot3d, [[cos(t)-1, y, -sin(t)], y = 1 .. 5, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

display([s0, ax, seq(s || i, i = 1 .. 3), s], scaling = constrained, lightmodel = light3, axes = normal, orientation = [60, 70], orientation = [70, 75], labels = [x, y, z]);

@nMaple 

restart;

with(plots):  with(plottools):

ax := spacecurve([-1, t, 0], t = -1 .. 6, color = red, thickness = 2, linestyle = dash):

s1 := animate(plot3d, [[(x+1)*cos(t)-1, x^2+1, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

s2 := animate(plot3d, [[(x+1)*cos(t)-1, 5, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

s3 := animate(plot3d, [[cos(t)-1, y, -sin(t)], y = 1 .. 5, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

display(ax, seq(s||i, i = 1 .. 3), scaling = constrained, lightmodel = light3, axes = normal, orientation = [60, 70], labels = [x, y, z]);

See my answer in this thread