Use plots:-shadebetween command.

Your example:
restart;
plots:-shadebetween(x+y+4, 25-x^2-y^2,x=0..1,y=x^2..2-x^2, showboundary = false, changefill = [color = khaki, transparency = 0], axes=normal, lightmodel=light1);

Addition. The simplest way to get the projection of a body onto some coordinate plane is simply to specify the corresponding orientation. For example, the projection of this body on the xOz plane (the x axis goes to the left, the z axis goes up):

 plots:-shadebetween(x+y+4, 25-x^2-y^2,x=0..1,y=x^2..2-x^2, showboundary = false, changefill = [color = khaki, transparency = 0], axes=normal, orientation=[90,90,0], lightmodel=none);

This animation clearly shows how the set of solutions (red segments) changes when the parameter a changes in the range  - 6 .. 0. 

A:=plots:-animate(plot,[x^3 + a*x + 2, x=-4..4, -15..15, color=blue, thickness=2], a=-6..0, frames=121):
B:=proc(a)
local sol, n, p1, p2;
uses plots;
sol:=fsolve(x^3 + a*x + 2);
n:=nops([sol]);
if n=3 then p1:=plot(0,x=sol[1]..sol[2], color=red, thickness=3);
p2:=plot(0,x=sol[3]..4, color=red, thickness=4) else
p2:=plot(0,x=[sol][1]..4, color=red, thickness=4) fi;
display(`if`(n=3,[p1,p2],p2), view=[-4..4,-15..15]);
end proc:
C:=plots:-animate(B, [a], a=-6..0, frames=121):
plots:-display(A, C);

Edit.

The first term is 1, and each next term is equal to the previous term written down by the same digits in the reverse order plus the position of the previous member.

The procedure for this:

P:=proc(n)
local L, N;
option remember;
if n=1 then return 1 else
L:=convert(thisproc(n-1), base, 10);
N:=nops(L);
n-1+add(L[i]*10^(N-i), i=1..N) fi;
end proc:

Example of use:

seq(P(n), n=1..100);
1, 2, 4, 7, 11, 16, 67, 83, 46, 73, 47, 85, 70, 20, 16, 76, 83, 55, 73, 56, 85, 79, 119, 934, 463, 389, 1009, 9028, 8237, 7357, 7567, 7688, 8899, 10021, 12035, 53056, 65071, 17093, 39109, 90232, 23249, 94273, 37291, 19316, 61435, 53461, 16481, 18508, 80629, 92657, 75679, 97708, 80831, 13861, 16885, 58916, 62041, 14083, 38099, 99142, 24259, 95303, 30421, 12466, 66485, 58531, 13651, 15698, 89719, 91867, 76889, 98938, 84061, 16121, 12235, 53296, 69311, 11473, 37489, 98552, 25669, 96733, 33851, 15916, 62035, 53111, 11221, 12298, 89309, 90487, 78499, 99578, 87691, 19771, 17885, 58966, 67081, 18173, 37279, 97372

Your system has infinitely many solutions, depending on one parameter (Maple has selected  y  as a parameter). Each equation of the system separately sets a plane in 3D. The  solve command easily finds all solutions as the intersection of these planes. In fact, we have a whole straight line of solutions (red line at the plot):

sol:=solve([x+y = 373320, z = (x+y) / 0.44 - y -  y* (1 - 0.99)], [x,y,z]);  # All the solutions
plots:-display(plots:-spacecurve(rhs~(op(sol)), y=0..400000, color=red, thickness=3),plots:-implicitplot3d([x+y = 373320, z = (x+y) / 0.44 - y -  y* (1 - 0.99)], x=0..500000,y=0..800000,z=0..1000000, style=surface,  color=[yellow,blue], axes=normal, tickmarks=[5,5,5]));  # Visualization

Use  Sum  instead of  sum:

f(x):=Sum((-1)^n*x^(2*n+1)/(2*n)!, n=0..infinity);
simplify(diff(f(x), x));

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.