Plot the anti derivative and remember the very theorem you want to use.

Int(sqrt(exp(2*I*t)-1),t);
J:=value(%);
[Re(%), Im(%)];
plot(%, t=0 .. Pi, color=[red, blue]);

You get it through the bivariate standard normal distribution N2 using
a transformation of variables:

restart;

f := proc (t1, t2, alpha1, beta1, alpha2, beta2, rho) options operator, arrow; 
  (1/4)*(sqrt(beta1/t1)+(beta1/t1)^(3/2))*(sqrt(beta2/t2)+(beta2/t2)^(3/2))*
  exp(-((sqrt(t1/beta1)-sqrt(beta1/t1))^2/alpha1^2+
  (sqrt(t2/beta2)-sqrt(beta2/t2))^2/alpha2^2-2*rho*(sqrt(t1/beta1)-
  sqrt(beta1/t1))*(sqrt(t2/beta2)-sqrt(beta2/t2))/(alpha1*alpha2))/(2-2*rho^2))/
  (alpha1*beta1*alpha2*beta2*Pi*sqrt(1-rho^2)) end proc;

N2:= (x,y,rho) -> 1/sqrt(1-rho^2)/2/Pi*
  Int(Int(exp(-(xi^2-2*rho*xi*eta+eta^2)/(2*(1-rho^2))), 
                    eta=-infinity..y), xi=-infinity..x);

trans := {x1 = (sqrt(t1/beta1)-sqrt(beta1/t1))/alpha1, 
          x2 = (sqrt(t2/beta2)-sqrt(beta2/t2))/alpha2};

assume(-1 < rho and rho < 1, 
       alpha1 > 0, beta1 > 0, alpha2 > 0, beta2 > 0, 
       t1 > 0, t2 > 0);

'N2(x1,x2,rho)': 'eval(%, trans)':
'diff(%, t1)': 'diff(%, t2)':;
'1/2*f(t1, t2, alpha1, beta1, alpha2, beta2, rho) = %';
 
simplify(%):
is(%);

  1/2 f(t1, t2, alpha1, beta1, alpha2, beta2, rho) =

           2
          d
        ------- eval(N2(x1, x2, rho), trans)
        dt2 dt1


                                 true

The task can be seen as to find Beta(1/n,n-1/n) = n, a fixed point. Using the integral theorem for the Beta function it means that the integrals over t = 0 ... 1 for  t^(1/n-1)*(1-t)^(n-1/n-1) and t^(1/n-1) are equal.

This certainly is the case if it holds true for the integrands. For that 'solve' includes the desired formal solution.

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