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));

@mweisbr 

So, you are using comma as a decimal separator. Probably this should work if the OS "locale" is set correctly.

But the simplest cure is to replace in your file "," by "." using a text editor (or Maple itself).

@Preben Alsholm 

in this case we still have:

Digits:=3:
[0.3333] <> [0.3333] + [0];

(at least visually.)

@Preben Alsholm 

The u in the post can be used. More generally:

# Proof
restart;
u:=evalf[Digits+1](1/3):
is( [u+0] = [u] + [0] ), evalb( [u+0] = [u] + [0] );

                            false, false

@Preben Alsholm 

There exists an algebraic exppresion u such that:

[u+0] <> [u]+[0];

@Cryme 

g is not C^2 (already answered).

A good package that does non need C^2 and finds even global optima is DirectSearch, see
https://www.maplesoft.com/applications/view.aspx?SID=101333

@mmcdara 

Because the question was about L^* . L = A.

@tomleslie 

But solve works directly. It is evalc which does not.

You have used evalc in your previous reply, and indeed

max(evalc(s));

works there.

@tomleslie 

Of course it is a symbolic problem. Maple should try harder to identify the real numbers, and not waiting for workarounds.
Note that evalc does not always work. Try:

s:=[solve(5*(x^2-1)*(x^2-x-4)+1,explicit)];
max(s);

The integration is not complicated. u is a polynomial/r^4, so the multiple integrals can be computed symbolically.
Use rationals instead of floats.
The problem is that you have many symbolic terms so the result will be huge.

@Robert Israel 

But the best definition is of course  :-)

f := n -> (n^2+n)^2/4;