I do not quite understand you first question, but for the second one:

sol:=RootOf(sin(x+y)+sin(x) = y, y);
taylor(sol, x = Pi, 15);

                             3                 5                   7
   - (x - Pi) + 1/12 (x - Pi)  - 1/240 (x - Pi)  + 1/10080 (x - Pi)

             17           9     1153            11      13297
         + ------ (x - Pi)  - --------- (x - Pi)   + -----------
           362880             159667200              24908083200

                13             15
        (x - Pi)   + O((x - Pi)  )

You can use simplify with side relations, see help for details on the simplify command

simplify(exp(a+b+c), {a+b+c=0});
                                  1

You have a typo, you need to write cos(...), not cos *

h:=t -> (2*t-1) * cos( sqrt(3*(2*t-1)^2+6) ) / (sqrt(3*(2*t-1)^2+6));

int(h(t), t);
                                 2            1/2
                    1/6 sin((12 t  - 12 t + 9)   )

PS: I converted your 'Post' to a 'Question'

You have several variables in your Z:

indets(Z, symbol);
                          {A0, A1, A2, r, x}

Of course you have to tell the system what is meant by Z', for example you can say
dZ:= diff(Z, x) or similar (instead of writing Z')

Re and Im depends on what you know about the variables. If you are sure they are
Reals then you may use evalc( Re(Z) ), same for Im. Else you just use Re and Im.

F:=eval(Integrand, [C = 1, beta = 1]);
MultiSeries:-series(F, s=0, 3);


shows, that it is like 1/s^2 in zero, an integral would be like 1/s, so it does not converge - the integral does not exist. Try starting in s=0.1 or similar.

Otherwise said: it is like asking "What to do for Int( 1/x^2, x = 0 .. infinity ) ?"

I think one can re-write it as

Sum( (r[j-1]*s+z)*exp(r[j]*z), j=1 .. 4) = - 4*z with r[0]:= r[4] where
d=s/L, z = beta z = beta*L and r = (indexed) Roots of x^4 + 4 = 0.

And suppose that it has infinitely many solutions for abs(z) --> infinity.

Strange ... I switched to exact values for S[j] and if I use w1 = 1/sqrt(1-x^2)
as weight and change to x = arcsin(t) then it plots pretty smooth (formally there
is a singularity in x=0), now over -Pi/2 .. Pi/2.

I would expect Maple to find it, at least after splitting in 0.

You could convert to a string and take its length

A:=x+2*y;
convert(A, string); length(%);
                                  5

To get rid of the length of variable names you can use LeafCount

B:=subs(x=theOneVariable, y=theOtherVariable, A);

               B := theOneVariable + 2 theOtherVariable

length(A), length(B);
                                9, 37

MmaTranslator[Mma][LeafCount](A), MmaTranslator[Mma][LeafCount](B);;

                                 5, 5

You can export through the menu bar File / Export as in various ways, HMTL, text formats (plain text or RTF) or similar through Save As (using 'Maple Input' generates a file *.mpl which is pure text) 

PS: Delete output within the sheet in case you want to see the input commands. That can be done through the menu bar Edit / Remove Output

There are two things.

In your integration over MI you mixed up the variables, try it step by step:

  H:= z -> evalf( Int( x -> MI(x,z), 0..3,method=_d01akc, epsilon=1e-8) );

  evalf( Int(H, 0 .. 2, method=_d01akc, epsilon=1e-8) );

                           1.43166543880904

Which differs a lot from the true value 24/Pi^2 ~ 2.413...

But the approximation is quite bad

  plot3d(sin(Pi*x/Lx)*sin(Pi*z/Lz) - 'MI'(x,z), 
    x=0..3,z=0..2, axes=boxed, title="absolute errors");

I used 'quotes' for MI to enforce simple plotting over a range which
is more convenient than feeding the approximation points.

PS: this is 2-dim integration, you integrate over 2 variables

tmp:=x^2+y^2 - 1:

indets(tmp, symbol); # ok, gives only w
indets(tmp);         # so you need 0 < w because of roots

plot(abs(tmp));      # w should be close to zero
plot(tmp, w=0 .. 1, 0 .. 10 );

Now try

fsolve(tmp, w=0 .. 1);
                       0.0240992450396199

Here is a very simple version for the bug: 1 / linear (please also see Carl's answer,
http://www.mapleprimes.com/questions/212935-Computing-A-Complex-Residue#answer228699).

So I guess the bug is not due to 'residue' itself, but some simplifications
or transformations.

 restart:
 g:=1/(z-sqrt(-1+I));
                                     1
                         g := ---------------
                                          1/2
                              z - (-1 + I)

 z0:=2^(1/4)*exp(3*Pi*I/8); 
 z1:= simplify( convert(%, radical) );

                             (1/4)
                      z0 := 2      exp(3/8 I Pi)

                    1/4        1/2 1/2           1/2 1/2
                   2    ((2 + 2   )    I + (2 - 2   )   )
             z1 := --------------------------------------
                                     2

 #eval(g, z =z1); radnormal(%); # div by zero -> z1 = pole
 z2:=sqrt(-1+I);

                                        1/2
                          z2 := (-1 + I)

 '0=z2 - z1'; evalc(%): simplify(%): is(%);
 '0=z2 - z0'; evalc(%): simplify(%): is(%);

                             0 = z2 - z1

                                 true

                             0 = z2 - z0

                                 true

 'residue(g, z = z0)': '%'=%;
 'residue(g, z = z1)': '%'=%;
 'residue(g, z = z2)': '%'=%;

                        residue(g, z = z0) = 0

                        residue(g, z = z1) = 0

                        residue(g, z = z2) = 1

 'residue(1/(z-a), z = a)': '%'=%;

                                1
                      residue(-----, z = a) = 1
                              z - a

 
PS: one can find that through looking at the partial fractions
and checking for which summand it fails

This is "by design" - consider the case that the integral can not be found.

You can use the following:

f:= x -> x;
F:= x -> int( f(X), X= 0 .. x);

F(x);
                                   2
                                  x
                                 ----
                                  2

# may be you want the following:

int( f(X), X= 0 .. x);
F1:=unapply(%, x);
                                   2
                                  x
                                 ----
                                  2

                                          2
                          F1 := x -> 1/2 x

"fsolve" gives a numerical result, depending on "Digits" - and it is a good idea to use Digits:=15 as default.

"solve" gives a symbolic result - but sometimes not all of them

You may use

solve( exp(z) = -1, allsolutions = true);

                          Pi I + 2 I Pi _Z1

Here the strange notation _Z1 is Maple's standard way to say that this is an arbitrary integer.

MP_compute_powers_loop.mws

Download MP_compute_powers_loop.mws

Maple has a "memory" for results, so times are not constant.

For computing it depends how you instruct Maple, for usual double precision that 1 million operations can be done from ~ 30 sec donw to 1/8 sec, see file attached