> is(0< mu<1);

(entered in 2D Math input) produces the result false  (i.e. Maple is just telling you that this does not follow from the assumptions you've made).  I think you want

> assume(0 < mu < 1);

But as for your error: it looks like you've assigned a value to theta(t)  in the line

theta(t):= ...

Once you've done that, you don't have theta(t) as an unknown function for your ODE any more.

Use textplot.

Use textplot.

So

2.*arctan(RootOf(3630444634*arctanh(5616816483/10000000000*_Z+8023282641/ 10000000000)-4012053071+50000000000*s0))+6.283185307*_Z1

means

2. * arctan(u) + 6.283185307*n

where n is an integer, and 3630444634*arctanh(5616816483/10000000000*u+8023282641/ 10000000000)-4012053071+50000000000*s0 = 0

So

2.*arctan(RootOf(3630444634*arctanh(5616816483/10000000000*_Z+8023282641/ 10000000000)-4012053071+50000000000*s0))+6.283185307*_Z1

means

2. * arctan(u) + 6.283185307*n

where n is an integer, and 3630444634*arctanh(5616816483/10000000000*u+8023282641/ 10000000000)-4012053071+50000000000*s0 = 0

Mathematical theories are not "registered", they are published.  If you have something that you think is worth publishing, you should find an appropriate mathematical journal.  Read the information for authors, which is generally on an inside cover or at the back of each issue.  What journal might be appropriate will depend on the subject matter and mathematical level of the material. 

The _Z and _Z1~ play completely different roles (and it was perhaps an unfortunate choice by Maple to use the same letter Z for both).  _Z in a RootOf is "the variable for the function that this is the root of", i.e. RootOf(_Z^2 - 2) means "some value z such that z^2 - 2 = 0".  There is no restriction on what type of complex number it could be.  On the other hand, _Z1~ means an arbitrary integer.

The _Z and _Z1~ play completely different roles (and it was perhaps an unfortunate choice by Maple to use the same letter Z for both).  _Z in a RootOf is "the variable for the function that this is the root of", i.e. RootOf(_Z^2 - 2) means "some value z such that z^2 - 2 = 0".  There is no restriction on what type of complex number it could be.  On the other hand, _Z1~ means an arbitrary integer.

In either case, I believe dsolve will use variable step sizes for solving the differential equation.  With the default settings, I don't think odeplot gives you adaptive plotting: in my experiments, the interval between the t values plotted was constant.  The only case where  odeplot gives you all the points actually calculated, I think, is if you give dsolve the range option and odeplot the refine option.  Thus

> sol:= dsolve({de, ic}, numeric,range=0..10);
  odeplot(sol, [t, y(t)], refine=1);

In either case, I believe dsolve will use variable step sizes for solving the differential equation.  With the default settings, I don't think odeplot gives you adaptive plotting: in my experiments, the interval between the t values plotted was constant.  The only case where  odeplot gives you all the points actually calculated, I think, is if you give dsolve the range option and odeplot the refine option.  Thus

> sol:= dsolve({de, ic}, numeric,range=0..10);
  odeplot(sol, [t, y(t)], refine=1);

I'm not sure where you're going with Schanuel's conjecture.  Lindemann's Theorem can tell you that the nonzero solutions to x = tan(x)  are transcendental: if x is algebraic and nonzero, so is ix, and then Lindemann says w = exp(ix) is transcendental.  But tan(x) = -i (w^2-1)/(w^2 + 1), so if tan(x) = x we'd have w = (+/-) sqrt((i-x)/(i+x)), which is algebraic if x is algebraic.

I'm not sure where you're going with Schanuel's conjecture.  Lindemann's Theorem can tell you that the nonzero solutions to x = tan(x)  are transcendental: if x is algebraic and nonzero, so is ix, and then Lindemann says w = exp(ix) is transcendental.  But tan(x) = -i (w^2-1)/(w^2 + 1), so if tan(x) = x we'd have w = (+/-) sqrt((i-x)/(i+x)), which is algebraic if x is algebraic.

Use odeplot in the plots package to plot the output of dsolve(..., numeric).  For example:

> exactsol:= dsolve({de, ic}, y(t));
   #  I assume that will produce an explicit formula y(t) = ... 
  numericsol:= dsolve({de, ic}, y(t), numeric);
  plots[odeplot](numericsol, [t, y(t) - rhs(exactsol)],t=0..50);

Use odeplot in the plots package to plot the output of dsolve(..., numeric).  For example:

> exactsol:= dsolve({de, ic}, y(t));
   #  I assume that will produce an explicit formula y(t) = ... 
  numericsol:= dsolve({de, ic}, y(t), numeric);
  plots[odeplot](numericsol, [t, y(t) - rhs(exactsol)],t=0..50);

> ListTools:-Flatten([[1,15,21],[5,7,33]]);

And then, as Georgios suggests, you can use sort on the result.