q(x, s) := -(-(-thetac*s^(3/2)-s^(5/2)+(s^2+s*thetac)*alpha1*sqrt(Dp))*Dc*A1*exp((lh-x)*sqrt(s+thetac)/sqrt(Dc))+((alpha1*(s+thetac)*(-pinf*s-pinf*thetac+A2)*Dp^(3/2)+s*sqrt(Dp)*(A1*(s+thetac)*sqrt(Dc)-Dc*alpha1*(-pinf*s-pinf*thetac+A2)))*sqrt(s+thetac)+A1*sqrt(Dp)*s*Dc*alpha1*(s+thetac))*exp(sqrt(s)*(-lh+x)/sqrt(Dp))-(-_F1(s)*(-s*alpha1*(s+thetac)^2*Dp^(3/2)-s^(3/2)*Dp*thetac^2+thetac*(Dc-2*Dp)*s^(5/2)+(Dc-Dp)*s^(7/2)+sqrt(Dp)*s^2*Dc*alpha1*(s+thetac))*exp((-2*lh+x)*sqrt(s)/sqrt(Dp))+_F1(s)*(-s*alpha1*(s+thetac)^2*Dp^(3/2)-thetac*(Dc-2*Dp)*s^(5/2)+(-Dc+Dp)*s^(7/2)+s^(3/2)*Dp*thetac^2+sqrt(Dp)*s^2*Dc*alpha1*(s+thetac))*exp(-sqrt(s)*x/sqrt(Dp))+alpha1*(s+thetac)*(-pinf*s-pinf*thetac+A2)*Dp^(3/2)+(-pinf*(Dc-2*Dp)*thetac+A2*(Dc-Dp))*s^(3/2)-pinf*(Dc-Dp)*s^(5/2)-s*alpha1*Dc*(-pinf*s-pinf*thetac+A2)*sqrt(Dp)-sqrt(s)*Dp*thetac*(-pinf*thetac+A2))*sqrt(s+thetac))/((s+thetac)^(3/2)*s*((Dc-Dp)*s-Dp*thetac)*(sqrt(Dp)*alpha1-sqrt(s))):

collect(q(x,s), [exp((lh-x)*sqrt(s+thetac)/sqrt(Dc)), exp(sqrt(s)*(-lh+x)/sqrt(Dp)), exp((-2*lh+x)*sqrt(s)/sqrt(Dp)), exp(-sqrt(s)*x/sqrt(Dp))]):

map(simplify, %);  # The final result

max[index]  gives only one solution. But the matrix can have multiple entries with the same maximum value. The following example shows how to find the positions of all of these entries:

restart;
with(LinearAlgebra):
p:=RandomMatrix(10,10, generator=0..12);
mval:=max[index](p):  # The position of one maximum entry of the matrix p
N:=p[mval];  # The maximum entry of the matrix p
k:=0:
for i to 10 do
for j to 10 do
if p[i,j]=N then k:=k+1; L[k]:=[i,j] fi;
od: od:
L:=convert(L,list);  # All the positions of maximum entries

For example:

r:=rand(-9..9):
k:=0:
do
a,b:=r(),r();
if igcd(a,b)=1 then k:=k+1; L[k]:=a*x+b fi;
if k=4 then break fi;
od:
convert(L, list);

In Help we can read "Note: The algsubs command currently works only with integer exponents". The positive integer exponents are here implied.

The toy example:
algsubs(x=t, x+1/x);
                                      t+1/x

In your the first example Maple simply makes constant  1/x  outside the integral and calculates the remaining integral. 

The second example is done correctly:

int(1/sin(x), x);
simplify(%);
convert(%, tan);
applyop(normal, 1, %);
                            

 or

int(1/sin(x), x);
diff(%, x);
simplify(%);
                                  

In the expression  u  in the last fraction in the numerator and in the denominator, sqrt(sigma)  is present. I removed it and now everything is working:

restart;
u := c-6*mu*(1+lambda*sqrt(-mu)*coth(A+sqrt(-mu)*(x+y+mu*t)))/((a+b)*sqrt(-mu)*coth(A+sqrt(-mu)*(x+y+mu*t)))+6*mu*sqrt((1-coth(A+sqrt(-mu)*(x+y+mu*t))^2))/((a+b)*sqrt(-mu)*coth(A+sqrt(-mu)*(x+y+mu*t)));

c := 1:                                     
mu := -1.5:                                     
lambda := 1.5:                                     
a := 1:                                      
b := 1:                                      
A := 0.5:                                     
y := 0:

plot3d(abs(u), x = -10 .. 3, t = -10 .. 3);

DeleteZeroLines:=proc(A::Matrix)
local m, n, ZR, L1, B, m1, ZC, L2;
uses LinearAlgebra;
m, n:=op(1, A);
ZR:=ZeroVector[row](n);
L1:=[seq(`if`(Equal(A[i], ZR), i, NULL), i=1..m)];
B:=DeleteRow(A, L1);
m1:=op([1,1], B);
ZC:=ZeroVector(m1);
L2:=[seq(`if`(Equal(B[..,i], ZC), i, NULL), i=1..n)];
DeleteColumn(B, L2);
end proc:

Example of use:

A:=RandomMatrix(8, 10, density = 0.1, generator = 1..9);
DeleteZeroLines(A);

                        

DeleteZeroLines.mw

Of course, this conversion can be done by direct calculation (without  convert  command):

z := arctan(x/sqrt(a^2-x^2)):
t:=op(1,z):
simplify(arcsin(t/sqrt(1+t^2)))  assuming a>0, x^2<a^2;

                                            arcsin(x/a)

convert((t-1)/(-1+2*t), parfrac);

                                 1/2-1/(2*(-1+2*t))

Replace the comma after  l2c2:=0  with a colon:

l2c1:= 2*k:  l2c2:=0:  l2c3:=8-k:
plots:-spacecurve([l2c1, l2c2, l2c3], k=-10..10, color=red);

A procedure for this:

restart;
MakeColored:=(s,c)->Typesetting:-mo(convert(s,string), mathcolor = c):

Examples of use:

x:=MakeColored(x, "Red"):
y:=MakeColored(y, "Blue"):
z:=MakeColored(z, "Green"):
x, y, z;
(x+y)/(x-z);

See these threads for details:

http://www.mapleprimes.com/questions/203931-Print-Hello-In-Green

http://www.mapleprimes.com/questions/126174-How-Do-You-Color-Text-Output-From-Command-Line

Edit.

 

I think that this is the case when it is easier to solve by hand, because it is easy to express exp(x) rationally through tanh(x/2). If you need to perform a similar conversion more than once, the procedure called  Convert  will be useful. It converts any occurrence of  exp  function into  tanh  function:

restart;
Convert:=proc(Expr)
local A;
convert(tanh(z), exp);
combine(normal(subs(exp(-z)=1/exp(z), %)));
solve(%=tanh(z), exp(2*z));
A:=subs(z=z/2, %);
simplify(applyrule(exp(z::anything)=A, Expr));
end proc:

Examples of use:

Convert(exp(x));
Expr:=-(exp(-a*s)-1)*kw/((exp(-a*s)+1)*s^2);
Convert(Expr);  # The initial example

Maple uses radians by default, so should be Pi/2 instead of 90.

Output := Array(-10 .. 10, -10 .. 10, (i, j)->i^2+j^2):
M:=convert(Output, listlist):
plots:-listcontplot(M, axes=none);

lambda:=2.6:
fsolve(sum(lambda^k*exp(-lambda)/k!, k=0..n)=0.95);

                                  4.980237188

Example:

restart;
sol:=solve({x-2*y=5, x+3*y=1});
eval(x, sol);  # The first way
eval(y, sol);
assign(sol);  # The second way
x, y;

If you need decimals, use  fsolve  instead of  solve.