The Maple dsolver is very powerful, but everything has advantages and disadvantages. I was recently asked the following question.
Let us consider the system of ODEs
>restart; sys := [diff(y(x), x) = -(4*cos(x)*y(x)+z(x)*cos(x)^2+3*z(x))/(sin(x)*(cos(x)^2-9)),
>diff(z(x), x) = -(y(x)*cos(x)^2+3*y(x)+4*z(x)*cos(x))/(sin(x)*(cos(x)^2-9))]:
The functions
>y1 := C[1]*(cos(x)+1)^(1/2)/(cos(x)+3)^(1/2)+C[2]*(1-cos(x))^(1/2)/(3-cos(x))^(1/2):
>z1 := -C[1]*(cos(x)+1)^(1/2)/(cos(x)+3)^(1/2)+C[2]*(1-cos(x))^(1/2)/(3-cos(x))^(1/2):
are its general solution in view of
>odetest({y(x) = y1, z(x) = z1}, [sys[1], sys[2]]);
  [0, 0]
However, the command
>sol := dsolve(sys, [y(x), z(x)]);
produces a long complicated output with the imaginary unit I. By the way, the  MTM[dsolve] command produces a similar output and Mathematica 9.0.0 simply returns the command DSolve[{y'[x] == -(4*Cos[x]*y[x] + z[x]*Cos[x]^2 + 3*z[x])/
(Sin[x]*(Cos[x]^2 - 9)), z'[x] == -(y[x]*Cos[x]^2 + 3*y[x] + 4*z[x]*Cos[x])/(Sin[x]*(Cos[x]^2 - 9))},
{y, z}, x].
The question is: how to obtain y1 and z1 in Maple?
Here is my answer. It is clear that all the simplifications of the output are a monkey business.
Having executed  the infolevel[dsolve]:=5:sol := dsolve(sys, [y(x), z(x)]); commands, I observed that place:
.......

 Change of variables used:
      [x = 1/2*arccos(t)]
   Linear ODE actually solved:
      (10*t-26)*u(t)+(478*t-6*t^3-70*t^2+110)*diff(u(t),t)+(-40*t+40*t^3-476+480*t^2-4*t^4)*diff(diff(u(t),t),t) = 0
........
That led me to believe that the reformulation of sys may be useful.
 Having remembered of the universal trigonometrical substitution, I tried
>simplify(eval(sys, [sin(x) = 2*t/(1+t^2), cos(x) = (1-t^2)/(1+t^2)]));

[diff(y(x), x) = -(1/2)*(1+t^2)*(-y(x)+y(x)*t^4-z(x)-z(x)*t^2-z(x)*t^4)/((2+5*t^2+2*t^4)*t),
diff(z(x), x) = (1/2)*(1+t^2)*(y(x)+y(x)*t^2+y(x)*t^4+z(x)-z(x)*t^4)/((2+5*t^2+2*t^4)*t)]

and after that
>eval(simplify(eval(sys, [sin(x) = 2*t/(1+t^2), cos(x) = (1-t^2)/(1+t^2)])), t = tan((1/2)*x));

[diff(y(x), x) = -(1/2)*(1+tan((1/2)*x)^2)*(-y(x)+y(x)*tan((1/2)*x)^4-z(x)-
z(x)*tan((1/2)*x)^2-z(x)*tan((1/2)*x)^4)/((2+5*tan((1/2)*x)^2+
2*tan((1/2)*x)^4)*tan((1/2)*x)), diff(z(x), x) = (1/2)*(1+tan((1/2)*x)^2)*(y(x)+
y(x)*tan((1/2)*x)^2+y(x)*tan((1/2)*x)^4+z(x)-z(x)*tan((1/2)*x)^4)/((2+5*tan((1/2)*x)^2+
2*tan((1/2)*x)^4)*tan((1/2)*x))]

>sol1 := dsolve(eval(simplify(eval(sys, [sin(x) = 2*t/(1+t^2), cos(x) = (1-t^2)/(1+t^2)])),
t = tan((1/2)*x))):
>sol1[1];

y(x) = (2*(_C1*sin((1/2)*x)*sqrt(6+2*cos(x))+
_C2*cos((1/2)*x)*sqrt(6-2*cos(x))))/(sqrt(6-2*cos(x))*sqrt(6+2*cos(x)))

Next steps were

>expand(rhs(sol[1]));

2*_C1*sin((1/2)*x)/sqrt(6-2*cos(x))+2*_C2*cos((1/2)*x)/sqrt(6+2*cos(x))

>simplify(applyrule([sin((1/2)*x) = ((1-cos(x))*(1/2))^(1/2), cos((1/2)*x) =
((cos(x)+1)*(1/2))^(1/2)], expand(rhs(sol[1]))), radical);

(_C1*sqrt(2-2*cos(x))*sqrt(6+2*cos(x))+
_C2*sqrt(2*cos(x)+2)*sqrt(6-2*cos(x)))/(sqrt(6-2*cos(x))*sqrt(6+2*cos(x)))

expand(%);

_C1*sqrt(2-2*cos(x))/sqrt(6-2*cos(x))+_C2*sqrt(2*cos(x)+2)/sqrt(6+2*cos(x))
 The happy end.
 I leave  the consideration of sol1[2] to the MaplePrimes users.

simplification.mw