g:=-18+2*x-8*y+5*z=0;
a:=tcoeff(lhs(g));
g-a;

                   -18 + 2 x - 8 y + 5 z = 0
                              -18
                      2 x - 8 y + 5 z = 18
 

U := {A[1], A[2], B[2]}:
op~(0, U);
                                   {A, B}

Addition. Another option with  map  command:

map(u->op(0,u), U);

AdvanceDate("Jan-01-1981", 13214);
AdvanceDate(%, -13214);

                                                     "Mar-07-2017"
                                                     "Jan-01-1981"

First, I plotted your equations in the desired ranges, they specified the boundaries for the roots:

Digits:=20:

plots:-implicitplot([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], beta=-0.1..1, delta=-0.1..1, color=[red,blue], thickness=2,  gridrefine=5, axes=box);
 
fsolve([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], {beta=0.2..0.6, delta=-0.1..0.1});  # The first solution 
fsolve([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], {beta=-0.1..0.1, delta=0.2..0.6});  # The second solution
 

 Addition. You can substitute the solutions found in the equations and make sure that they satisfy them with high accuracy.

In fact, the result is correct, but not simplified for x>0. I do not know a simple solution to the problem. To simplify the result to the desired form, I was able only by changing the variable:

restart;
A:=Int(1/(x+2*sqrt(x)), x);  # The initial integral in inert form
IntegrationTools:-Change(A, x=t^2, t);  # The change of variable x
simplify(combine(value(%), ln, anything, symbolic)) assuming t>0;  # The simplification
subs(t=sqrt(x), %);  # The reverse change 

restart;
int(2-8*x^2+32*x^4-128*x^6+512*x^8-2048*x^10, x);
sort(%, x, ascending);
 

The domain of the function  z->arctanh(z)  in real domain is  -1<z<1 . See

FunctionAdvisor(definition, arctanh);

But  (exp(2*y)+sqrt((exp(2*y))^2+exp(2*y)))/exp(2*y)-1>1  for any real y .
 

Your system has the parameter  V__system. To obtain a numerical solution, you must specify a numerical value for this parameter. I took  V__system=100 :

sys := {I__Arc = V__system/sqrt(2*Z1+R__Arc), R__Arc = (8750+C)/I__Arc^1.4}:
C := 300:  Z1 := 4:  V__system := 100:
solve(sys, {I__Arc, R__Arc});

                       {I__Arc = 1.178819336, R__Arc = 7188.237669}

fsolve  without additional options fails with this system.

For an expansion of the initial expression  ee  I slightly shortened acer's code and used expand  command with needed options:

ee:=sqrt(Dp)*(-Dp*sqrt(s+thetac)*alpha1*pinf*s^2-2*Dp*sqrt(s+thetac)*alpha1*pinf*s*thetac-Dp*sqrt(s+thetac)*alpha1*pinf*thetac^2+A2*Dp*sqrt(s+thetac)*alpha1*s+A2*Dp*sqrt(s+thetac)*alpha1*thetac+Dc*sqrt(s+thetac)*alpha1*pinf*s^2+Dc*sqrt(s+thetac)*alpha1*pinf*s*thetac+A1*Dc*alpha1*s^2+A1*Dc*alpha1*s*thetac+A1*sqrt(Dc)*sqrt(s+thetac)*s^2+A1*sqrt(Dc)*sqrt(s+thetac)*s*thetac-A2*Dc*sqrt(s+thetac)*alpha1*s)*exp((-lh+x)*sqrt(s)/sqrt(Dp))/((s+thetac)^(3/2)*s*(Dc*s-Dp*s-Dp*thetac)*(-sqrt(Dp)*alpha1+sqrt(s))):

new := simplify(radnormal(ee), size);

expand(new, exp((-lh+x)*sqrt(s)/sqrt(Dp)), s+thetac, pinf*s+pinf*thetac-A2, (Dc-Dp)*s-Dp*thetac);  # The final result 

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