e2_3:= w^(1-sigma)*((f__11*sigma*Delta__1/L__1)^(-1/(sigma-1))/w)^(k-sigma+1)*k/(Delta__1*(k-sigma+1)*a__0^k);
simplify(e2_3, symbolic);

hf:=<1/(1-x1-x2)^2+1/x1^2, 1/(1-x1-x2)^2; 1/(1-x1-x2)^2,1/(1-x1-x2)^2+1/x2^2>; # Matrix
x:=<10, 20>; # Vector
subs([x1=x[1], x2=x[2]], hf);

You wrote  x:=[10, 20];  This  x  is a list not a vector
 

plot([Re,Im]~([fsolve(x^3 = 1, complex)]), style=point, color=red, symbol=solidcircle, symbolsize=20, scaling=constrained);

Explanation. In fact, in one line we combined 3 steps.

The first step is to find three roots:

L:=[fsolve(x^3 = 1, complex)];
L := [-.500000000000000-.866025403784439*I, -.500000000000000+.866025403784439*I, 1.]
 

The second step is to extract the real and imaginary parts from each root:

[Re,Im]~(L);
[[-.500000000000000, -.866025403784439], [-.500000000000000, .866025403784439], [1., 0.]]
 

The last step is the plotting of these three points by plot command.

First we solve the system if in the first equation we replace the inequality by equality. The second equation is squared (we add the condition  -2*x^2+y^2>=0  for preserving the equivalence):

solve({y = 4*x^4+4*x^2*y+1/2, (1/2)*(x-y)^2-(x-y)^4 = (-2*x^2+y^2)^2, -2*x^2+y^2>=0}, [x, y]);

                                        [[x = -1, y = -3/2], [x = 0, y = 1/2]]

As for the initial system, all solutions will be obtained if we add to the points obtained all the points of the blue curve lying in the green region:

A:=plots:-implicitplot([y = 4*x^4+4*x^2*y+1/2, (1/2)*(x-y)^2-(x-y)^4 = (-2*x^2+y^2)^2, -2*x^2+y^2=0], x=-3..3, y=-3..1, color=[red,blue, green], thickness=2, gridrefine=5):
B:=plots:-inequal({-2*x^2+y^2>=0, y<=0},x=-3..3, y=-3..1, color="LightGreen"):
plots:-display(A, B, scaling=constrained);  

What do you want to get as

is an erroneous result, which is easily verified by differentiation. On the other hand, Maple calculates everything correctly.

See  

Download ccc_new.mw

The solution in two steps:

restart;
with(IntegrationTools):
M:=Int(diff(B(x,y),x,x)*F(x,y),[x,y,z])+Int(diff(E(x,y),x,x)*H(x,y),[x,y,z])+Int(diff(E(x,y),y,y)*F(x,y),[x,y,z])+Int(diff(B(x,y),x,x)*H(x,y),[x,y,z])+Int(diff(F(x,y),x,x)*F(x,y),[x,y,z])+Int(diff(E(x,y),x)*F(x,y),[y,z])+Int(diff(E(x,y),y)*H(x,y),[y,z])+Int(diff(B(x,y),y)*F(x,y),[y,z]):
M1:=map(t->applyop(collect,1,t, [F(x,y),H(x,y)]), Combine(M));
map(Int,op([1,1],M1),[x,y,z])+map(Int,op([2,1],M1),[y,z]);  # The final result 

The code was edited.

If you want to get a relationship between  x  and  a, b, p and r, then you just need to exclude the variables  y  and  z  from this system:

restart;
sys:={y=p*x, x=100*z/(p*r+r), z=x*a+y*(1-b)};
eliminate(sys,{y,z});
              [{y = p*x, z = (1/100)*x*r*p+(1/100)*x*r}, {x*(-100*b*p-p*r+100*a+100*p-r)}]

Thus, we have relation

x*(-100*b*p-p*r+100*a+100*p-r)=0

That is, x = 0 in the general case or if  -100*b*p-p*r+100*a+100*p-r=0  then  x  can be  an arbitrary number. 

I know only one way (though I never use it). After you have finished the typing of your code (the cursor must remain within the execution group), just click the right mouse button and select "convert to 2d math"

Example.
Here is the code of a simple procedure that searches for the N first pairs of prime twins.

The initial code:

restart;
PrimeTwins:=proc(N)
local n, a, b, L;
n:=0: a:=2:
do
b:=nextprime(a);
if b-a<=2 then n:=n+1; L[n]:=[a,b]; a:=b else a:=b fi;
if n=N then break fi;
od;
convert(L,list);
end:

After conversion:

a:=[1,2,3,4,5];
select(`<`,a,4);

restart;
with(IntegrationTools):
A:=diff(f(x,y),x)*diff(c(x,y),x);
B:=Int(Parts(Int(A,x), op(1,A)), [y,z]);
C:=Expand(B);

g:=(u,v)->alpha(u,v)*beta(u,v);
h:=(u,v)->alpha(u,v)+beta(u,v);
f:=(u,v)->D[1](g)(u,v)*h(u,v);
f(u,v);
                          

Or at once:

f:=unapply(D[1](g)(u,v)*h(u,v), u,v);      

Your system is an overdetermined one, because it has 2 unknown functions, but 3 equations. I just removed the last equation and now everything is all right. You can solve the last equation separately, because it depends only on g(eta).

sol1 := dsolve([diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2 = 0, diff(diff(diff(g(eta), eta), eta), eta)+f(eta)*(diff(diff(g(eta), eta), eta))+2*g(eta)*(diff(diff(f(eta), eta), eta))-3*(diff(g(eta), eta))*(diff(f(eta), eta)) = 0, f(0) = 1, (D(f))(0) = 1, (D(f))(5) = 0, g(0) = 1, (D(g))(0) = -1, (D(g))(5) = 0], numeric, method = bvp);
plots[odeplot](sol1, [[eta, f(eta)], [eta, g(eta)]], eta = 0 .. 5, color = [red, blue], scaling = constrained);

To see the solution steps, you can use this command:

Student:-Basics:-LinearSolveSteps( "-(2*x+10)+8=4*(-x+1)", x);

Use  solve  instead of  fsolve :

eqn1:=-x*y^2+4*x=5;
eqn2:=(1/3)*x^3+y^2=1;

sol:=solve({eqn1, eqn2, x>=-10, x<=10, y>=-10, y<=10}, {x,y});
evalf(%);

For plotting use  plots:-implicitplot  command.

f := L->max(abs~(L));

Example of use:

f([1+I, 2-I, -1+3*I]);
                                           10^(1/2)