plots:-implicitplot(w=-4.918231636*z^3+2.486861786*z^2+.8573483099*z+2.341623514, w=-3..7, z=-3..3);

I did not quite understand the question 2). If you just want a part of the graph from above in the specified ranges, then write:

plots:-implicitplot(w=-4.918231636*z^3+2.486861786*z^2+.8573483099*z+2.341623514, w=-2..3, z=-0.5..0.5, view=[-2..3, -0.5..0.5]);

If you want the initial equation (w against z) in the specified ranges, then write:

plot(-4.918231636*z^3+2.486861786*z^2+.8573483099*z+2.341623514, z=-2..3, -0.5..0.5);

if you have a new variable  u=z/H  then the change  z=H*u  should be done first with a specific  H  in the initial equation  w=-4.918231636*z^3+2.486861786*z^2+.8573483099*z+2.341623514

Edit.

In the help on  evalc  command we can read "The fundamental assumption that evalc makes is that unknown variables represent real-valued quantities". Therefore evalc immediately treats  z  as a real number and often ignores  assuming  or  assume .  When using  evalc with symbolic complex numbers just write  x+I*y  instead  z :

z:=x+I*y:
evalc(Re((Re(z)+I*Im(z))^2));
evalc(z-Re(z)-I*Im(z)); 

                                                x^2-y^2
                                                     0

Addition. If you want the result to be expressed through z-variable, you can do the reverse substitution:

z:=x+I*y:
evalc(Re((Re(z)+I*Im(z))^2)):
subs([x=Re('z'),y=Im('z')], %);

                                   Re(z)^2-Im(z)^2

Use the straightforward and simple syntax for this and you will not have any problems:

solve(0.042^x/0.021^x=4, x);
solve(0.260^y/0.065^y=1, y);

                          2.000000000
                               0.
 

In Maple 12 there was not  plottools:-getdata  command. Use  op  command instead. Your plot (let's call it P) contains several graphs.  op([1,1], P)  returns the data for the first graph as a numeric matrix, op([1,2], P)  is for the second graph and so on. Here is a simple example:

P:=plot([sin(x), cos(x)], x=-Pi..Pi, color=[red,blue]);
Q:=op(P);
data:=op([1,1], P);
plot(data);

Addition. If necessary, you can easily convert the received data from the matrix form to the list of lists, in which each sublist is the coordinates of some point on the graph:

convert(data, listlist);

Example:

restart;
A:=<a,b; c,d>;
LinearAlgebra:-Diagonal(A);  # The output is a column vector
# or
seq(A[i,i], i=1..op([1,1], A));  # The output is a sequence

Do not use linalg package. It's deprecated.

Also you can use  lhs  and  rhs  commands for this:

CI:=OneSampleZTest(S, mu, sdev, confidence=0.95, output='confidenceinterval');

CI_lower:=lhs(CI);

CI_upper:=rhs(CI);

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