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)

evalf~([solve(z^5-5*z^4-z^2-2)]);
[argument, abs]~(%);

But it is a correct answer:

limit(n^(3/2)*sum(2*k/(2*k+3), k=1..n), n=infinity);

                                    infinity

In Maple  e  is just a symbol, not  2.71828...

Should be:

int(exp(-x)*sin(x^2)/(2+x), x=1..infinity);
evalf(%);

Also, do not use square brackets to group expressions. They serve to form lists.

If you already know the result, you can just do :

is(abs(u+v)^2 + abs(u-v)^2 = 2*abs(u)^2 + 2*abs(v)^2);
                                          true

Your matrix equation can easily be reduced to a linear system:

gc();
restart:
with(LinearAlgebra):
A := <<a__11|a__12>,<a__21|a__22>>;
P := <<p__11|p__12>,<p__12|p__22>>;
Id := <<1|0>,<0|1>>;
eqn := Transpose(A).P+P.A =~ -Id;
solve({seq(seq(eqn[i,j], j=1..2), i=1..2)}, {p__11,p__12,p__12,p__22});
 

Although your system is linear with respect to the unknowns, it has 10 parameters. Even if Maple wrote formal expressions for the roots through the determinants, we would get huge expressions useless for use. I remind you that the matrix determinant 15 by 15 has in general 15! = 1307674368000 terms. Therefore, it is reasonable to solve the system by specifying the values of the parameters. See the example in the file.

System_new.mw