Download yourU.mw

@vv 

OK, I was looking only for symmetric solutions and I see now that your U is not symmetric.
Changing

U:=< u11,u12,0,  0;     # look for a block-diagonal U, such as M and C
     u21,u22,0,  0;
     0,  0,  u33,u34;
     0,  0,  u43,u44>;

==> 64 solutions, some of them depending on a parameter.
Most probably your U is one of them. Just check.

@Adam Ledger 

It is faster, and maybe e.g. the user only wants to have the duplicates near each other.

@Kitonum 

A not so obvious fact is that the sequence X(n) is dense in the interval [1,2], i.e. for each 1 <= t <= 2  there is a subsequence converging to t.

@das1404 

The lines can be removed. E,g.
display(polygon(A), axes=none, color=gold, style=patchnogrid);

@acer 

The real option in solve does not seem to be reliable.

restart;
solve((x^2-1)*(exp(x)+1),{x});
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},real);
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},real,allsolutions);
                 {x = I Pi}, {x = 1}, {x = -1}
solve((x^2-1)*(exp(x)+1),{x},allsolutions);
           {x = I Pi + 2 I Pi _Z1}, {x = 1}, {x = -1}

identify is far for being reliable. It fails here for some polynomials.

@Mariusz Iwaniuk 

@nyarko 

After all, sqrt is the simplest function after polynomials.
Just think: can the sqrt be removed in  2 + sqrt(3) ?

The method was posted just for fun and to show a (very inefficient!) possibility. The polygonal approach is much more efficient but designing nice letter shapes is time consuming. The ideal solution would be to have an angle parameter in textplot.

Or, including the asymptotics:

with(IntegrationTools):
A:=Int(exp(-x*t)/sqrt(t^2+t),t=1..infinity):
to 6 do
  A:=Parts(A,GetIntegrand(A)/exp(-x*t)) assuming x>0;
od:
eval(A,Int=0);

@rlopez 

@Math Pi Euler 

Probably you have included the prompt.
Here is the worksheet: asy.mw

Forget about the robot, focus on maths.

So, what are the equations (constraints), and what do you want to maximize?

If you only have 2 equations
cos(q[1])*cos(q[3])-sin(q[1])*sin(q[3])+sin(q[1])*q[2]+cos(q[1]) - 1 = 0,
sin(q[1])*cos(q[3])+cos(q[1])*sin(q[3])-cos(q[1])*q[2]+sin(q[1]) -1 = 0;

then there are infinitely many solutions. One of them is:

fsolve([
cos(q[1])*cos(q[3])-sin(q[1])*sin(q[3])+sin(q[1])*q[2]+cos(q[1]) - 1 = 0,
sin(q[1])*cos(q[3])+cos(q[1])*sin(q[3])-cos(q[1])*q[2]+sin(q[1]) -1 = 0,
q[1]=0
]);

        {q[1] = 0., q[2] = 0., q[3] = 1.570796327}

@Adam Ledger 

See the help page  ?semilogplot

Actually you can simply use  the equivalent 

plot(eval(u, x=10^t), t=1..100);

the only difference being the tickmarks.

A very nice post!
Unfortunately it shows that to obtain the desired asymptotic expansion
the user must have sound knowledge in several mathematical fields such as
Special functions, Complex analysis, Optimization theory etc.
And of course Maple programming.
Let's hope that in the near future a simple call to
series(), MultiSeries:-series() or maybe MultiSeries:-multiseries
will be enough in most cases.

Just a side note. Some math parts of the worksheet should be converted into
"non-executable" for easier step-by-step execution.

The type of evalf[j](hf)  is hfloat  for 10 <= j <= 15.

restart;
hf:=HFloat(1.234567890123456789);
type(hf,hfloat);
for j to 20 do
lprint(j,evalf[j](hf));
od:

                        hf := 1.23456789012346
                              true

1, 1.
2, 1.2
3, 1.23
4, 1.235
5, 1.2346
6, 1.23457
7, 1.234568
8, 1.2345679
9, 1.23456789
10, HFloat(1.23456789012345669)
11, HFloat(1.23456789012345669)
12, HFloat(1.23456789012345669)
13, HFloat(1.23456789012345669)
14, HFloat(1.23456789012345669)
15, HFloat(1.23456789012345669)
16, 1.234567890123457
17, 1.2345678901234567
18, 1.23456789012345669
19, 1.23456789012345669
20, 1.23456789012345669

So, a workaround for a correct evalf[11](hf)  would be e.g.
evalf[11](evalf[20](hf));