If the system contains parameters, the  LinearSolve  command does not give the correct result, even if all the dimensions are specified.

Example:

A:=Matrix([[a,1], [1,b]]):

B:=Vector([1,2]):

LinearAlgebra[LinearSolve](A, B);

For example, if  a=1/2  and  b=2, the result is incorrect.

with(plots):

D124 := implicitplot3d(max(-x+y+z, x-y+z, x+y-z) = 1, x = 0 .. 1, y = 0 .. 1, z = 0 .. 1, style=surface, numpoints = 50000):

D134 := implicitplot3d(max(-x-y-z, -x+y+z, x-y+z) = 1, x = -1 .. 0, y = -1 .. 0, z = 0 .. 1, style=surface, numpoints = 50000):

D123 := implicitplot3d(max(-x-y-z, -x+y+z, x+y-z) = 1, x = -1 .. 0, y = 0 .. 1, z = -1 .. 0, style=surface, numpoints = 50000):

D234 := implicitplot3d(max(-x-y-z, x-y+z, x+y-z) = 1, x = 0 .. 1, y = -1 .. 0, z = -1 .. 0, style=surface, numpoints = 50000):

display(D124, D134, D123, D234, axes=normal, view=[-1.7..1.9, -1.7..1.9, -1.7..1.9], orientation=[30, 50], lightmodel=light4);

When you write  contours=n  then maple  finds the corresponding values ​​of the function by the formulas

seq(m+(1/(n+1))*(M-m)*k,  k = 1 .. n) , 

in which  m  and  M  are the maximum and minimum values of the  function in the specified region.

Your example:

with(plots):

z := x^5-y^3:

m := minimize(z, x = -5 .. 5, y = -4 .. 3); M := maximize(z, x = -5 .. 5, y = -4 .. 3);

L := [seq(evalf[4](m+(1/11)*(M-m)*k), k = 1 .. 10)];

A := textplot([seq([surd(L[k]+(-3+.5*k)^3, 5), -3+.5*k, L[k], align = right], k = 1 .. 10)]):

B := contourplot(z, x = -5 .. 5, y = -4 .. 3, axes = boxed, filledregions = true, coloring = [yellow, red], contours = 10):

display(A, B);

A nonlinear equation can have any number of zeros from 0 to infinity. The  evalf  command finds only one root
Therefore,  firstly it is  useful to examine the properties of the function  s->tanh(s)-C3*s. Here are
two basic properties:
 
1) It is odd.
2) On the positive s-axis the function is strictly convex upward.

From these properties it follows that if  0 <C3 <1, then there exists a unique positive root  s> 0. Altogether there will be three roots  -s, 0, s .

The solution of this equation can be issued in the form of a procedure in which the formal arguments are  C3  and  N (number of digits in the result).

Sol := proc (C3, N)

 local sol, h, A, B;

Digits := N;

if 1 <= C3 or C3 <= 0 then print(`One solution`);  print(s = 0);

plot([tanh(s), C3*s], s = -1.5 .. 1.5, -1.5 .. 1.5, color = [red, blue]) else sol := fsolve(tanh(s) = C3*s, s = 0 .. infinity);

print(`Three solutions`); print(s[1] = -sol, s[2] = 0, s[3] = sol);

h := eval(tanh(s), s = sol);

A := plots:-textplot([[-sol, 0, s[1]], [0, 0, s[2]], [sol, 0, s[3]]], align = [left, above]);

B := plot([tanh(s), C3*s, [-sol, t, t = -h .. 0], [sol, t, t = 0 .. h]], s = -sol-1 .. sol+1, thickness = [2, 2, 1, 1], color = [red, blue, black, black], linestyle = [1, 1, 2, 2], `if`(sol < 5, scaling = constrained, NULL));

plots:-display(A, B);

end if;

end proc;

Example:

Sol(0.18, 20);

f := x^5+a*x^4+b*x^3+c*x^2+d*x+e;
algsubs(x^3=1, f);

                   f := x^5+a*x^4+b*x^3+c*x^2+d*x+e
                            (1+c)*x^2+(a+d)*x+b+e

Procedure called  P  is plotting  your function so that x-axis to be plotted  at y=y0 . Parameters  a  and  b  specify the range on x-axis.

P := proc (expr, a, b, y0)

local m, M;

m := min(minimize(expr-y0, x = a .. b), -1);

M := max(maximize(expr-y0, x = a .. b), 1);

plot(expr-y0, x = a .. b, y = m .. M, tickmarks = [default, [seq(i = i+y0, i = floor(m) .. ceil(M))]]);

end proc:

Example:

P(x^2-4, -2, 2, -5);

Edited. The x-axis can be anywhere: below, above or to cross the plot of the function. 

I copied your expression without the last term, as it is trimmed in your message.

A:=-4*cos((1/5)*Pi)*m^3-4*cos((1/5)*Pi)*n^3-4*cos((1/5)*Pi)*p^3-4*cos((1/5)*Pi)*q^3+4*cos((2/5)*Pi)*m^3+4*cos((2/5)*Pi)*n^3+4*cos((2/5)*Pi)*p^3+4*cos((2/5)*Pi)*q^3-10*cos((1/5)*Pi)*p*q^2+10*cos((2/5)*Pi)*m^2*n+12*cos((2/5)*Pi)*m^2*p+12*cos((2/5)*Pi)*m^2*q+12*cos((2/5)*Pi)*m*n^2+10*cos((2/5)*Pi)*m*p^2+12*cos((2/5)*Pi)*m*q^2+12*cos((2/5)*Pi)*n^2*p+10*cos((2/5)*Pi)*n^2*q+12*cos((2/5)*Pi)*n*p^2+12*cos((2/5)*Pi)*n*q^2+12*cos((2/5)*Pi)*p^2*q+10*cos((2/5)*Pi)*p*q^2-10*cos((1/5)*Pi)*m^2*n-12*cos((1/5)*Pi)*m^2*p-12*cos((1/5)*Pi)*m^2*q-12*cos((1/5)*Pi)*m*n^2-10*cos((1/5)*Pi)*m*p^2-12*cos((1/5)*Pi)*m*q^2-12*cos((1/5)*Pi)*n^2*p-10*cos((1/5)*Pi)*n^2*q-12*cos((1/5)*Pi)*n*p^2-12*cos((1/5)*Pi)*n*q^2-12*cos((1/5)*Pi)*p^2*q-24*cos((1/5)*Pi)*m*n*p-24*cos((1/5)*Pi)*m*n*q-24*cos((1/5)*Pi)*m*p*q-24*cos((1/5)*Pi)*n*p*q+24*cos((2/5)*Pi)*m*n*p+24*cos((2/5)*Pi)*m*n*q+24*cos((2/5)*Pi)*m*p*q:

simplify(convert(A, radical));

Example:

M := plottools:-arc([0, 0], 1, (1/6)*Pi .. 3*Pi*(1/4), color = blue, thickness = 4):

plots:-display(plot(op(M)[1], color = pink, filled, scaling = constrained), M);

Your sequence can be defined by recurrent formula or  procedure-function.

By recurrent formula:

By procedure:

a:=proc(n)

local b, i;

b[1]:=1; b[2]:=2;

for i from 3 to n do

if type(i, odd) then b[i]:=b[i-2]*b[i-1] else b[i]:=b[i-3]+b[i-2] fi;

od;

b[n];

end proc: 

Example:

seq(a(n), n=1..20);

1, 2, 2, 3, 6, 5, 30, 11, 330, 41, 13530, 371, 5019630, 13901, 69777876630, 5033531, 351229105131280530, 69782910161, 24509789089304573335878465330, 351229174914190691

restart;

with(plots):

A:=contourplot(x^3-y^2,x=-3..3,y=-3..2, thickness=2, axes=boxed, contours=[seq(-4..5)]):

B:=textplot([seq([surd(0.7^2+i,3),0.7, i], i=-4..5)], font=[TIMES,ITALIC,12], color=brown, align={left, above}):

display(A,B);

There is no  labelledcontourplot  command in Maple.

plots[contourplot](x^3-y^2,x=-3..3,y=-3..2, thickness=2, axes=boxed, contours=10);

3 options are working:

L:=[sum, value, plot]:

for i in L do

Student[Calculus1][ApproximateInt](sin(Pi*x), x=0..1, functionoptions=[filled=true, color=pink], output=i);

od;

If "to reduce"  means reduction in the number of lines of code, then possible variant

simplify(solve(map(t->t<2, [solve(x^4 -(3*m+2)*x^2 + 3*m+3=0, x)]), m));

                              {2/3*2^(1/2) <= m, m < 11/9}

To improve the quality of plot use  numpoints  option, for example

DEplot({de1, de2}, [x(t), y(t)], t = 0 .. 5, [[x(0) = 10, y(0) = 0]], stepsize = 0.5,

color = black, linecolor = blue, numpoints=5000);

For the second and third plots arrows are not drawn because the equations are non-autonomous.

Since solutions of  your system  [Er1, Ez1]  are symmetric respect to the line  rho1=0 , then we assume  rho1>0 . It is easy to check that the second equation has a solution   {z1=R/2, rho1 - any number} ,  that will be unique for  R>0  and approximately  R<=1.2 .For such values ​​of  R  the second unknown  rho1  can be obtained from the solution of the first equation for  z1=R/2 . This solution will be called trivial. The blue line on the Carl Love's plot  in the range  0<R<=1.2  is exactly a straight line with equation z1=R/2 .

If  R> 1.2, then except the trivial solution there are two non-trivial solutions.  Carl Love's code finds smaller solution  z1[1]<R/2 . Since  non-trivial solutions are symmetric respect to the line  z1=R/2 , then the second non-trivial solution be  z1[2]=R - z1[1] . 

All of the above conclusions are clearly seen from the plots of the equations of the system. Plots are made for specific  R  in the range  0.2<=R<=2  with step=0.2 . 

Branching_of_solutio.mw