@David Sycamore  You do not need to do anything with the code of this procedure. You just have to copy this code as plain text into your worksheet. Then you press the enter key and after that you can use the procedure as shown in the examples above. You type (or copy) the procedure name and type the desired prime number in parentheses.

@JAMET   I did not understand the meaning of the question.

@Carl Love   for this useful clarification.

@Carl Love  I know. Of course, if for OP is necessary to select the smallest elements, the sum of which does not exceed a bound, then we can use a set. Otherwise, we should use a list.

@Rouben Rostamian   If you take  n  large enough, (n-5)/6^6  will be greater than 1.

1. 4 is the the total number of the roots of the equation.

2. To solve this equation it is better to use  Student:-Calculus1:-Roots  command. In your example, it immediately returns all 4 roots in symbolic form:

Student:-Calculus1:-Roots((x-1)*(x^3-9*x^2+4), x);

@weidade37211  As a part of your equation, write  (1-10^(-15))^(10^10-1);  in a separate line and press enter (you get  overflow). This is a rational number with a giant numerator and denominator. Maple tries to calculate it exactly and overflow occurs. To avoid this, it is enough to put a decimal point on one of the exact numbers. Then all calculations will occur in float format with the accuracy established by Digits (with 10 digits by default).

(1-10^(-15))^(10^10-1);
Error, numeric exception: overflow

(1-10^(-15))^(10^10.-1);
                                                   1.

To avoid loss of accuracy, I increased  Digits  in my answer below. See:
Digits:=20:
(1-10^(-15))^(10^10.-1);
                                                 0.99999000005000083332

@Wavish  Since you do not have a clear criterion for obtaining the wished line, then the solution "by eye" is likely to be the easiest and best solution to the problem.

An example of the creation of a specific matrix function:

F:=unapply(<t, t^2; sin(t), cos(t)>, t):
F(2);

@vv  In fact, as a formal symbolic result, we get the same thing. But if we try to assign some numeric constants to the symbols, then of course we get an error:

Ec := (Ems+I*Eml)*(1+((Ems+I*Eml)/Ef-1)*Zeta*phi/((Ems+I*Eml)/Ef+Zeta))/(1-((Ems+I*Eml)/Ef-1)*phi/((Ems+I*Eml)/Ef+Zeta)):
a:=simplify(Re(Ec)) assuming positive;
b:=simplify(Im(Ec)) assuming positive;
EC:=eval(Ec,Zeta=ZETA):
simplify(evalc([Re,Im](EC)));
eval(%, ZETA=Zeta);
is(a=%[1] and b=%[2]);
Zeta:=1;

@bsoudmand  This is the output of my code (in Maple 2018.2):

In the text form:

a:=(-Ems*Zeta*(Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+((Zeta-1)*Ems^3+Ef*(Zeta-1)^2*Ems^2-(Zeta-1)*(Ef^2*Zeta-Eml^2)*Ems-Ef*Eml^2*(Zeta+1)^2)*phi+Ems*(Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2))/((Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(-2*Ems^2-2*Ef*(Zeta-1)*Ems+2*Ef^2*Zeta-2*Eml^2)*phi+Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2);

b:=-Eml*(Zeta*(Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(Ef*(Ef-2*Ems)*Zeta^2+(-Ef^2-Eml^2-Ems^2)*Zeta-2*Ef*Ems+Eml^2+Ems^2)*phi-Ef^2*Zeta^2-2*Ef*Ems*Zeta-Eml^2-Ems^2)/((Ef^2-2*Ef*Ems+Eml^2+Ems^2)*phi^2+(2*Ef*(-Ems+Ef)*Zeta+2*Ef*Ems-2*Eml^2-2*Ems^2)*phi+Ef^2*Zeta^2+2*Ef*Ems*Zeta+Eml^2+Ems^2);

@vv 

SymFun(ln((x-1)/(x+1)), x);
plot(ln((x-1)/(x+1)), x=-4..4, -4..4);

SymFun(sqrt(x^2), x);
plot(sqrt(x^2), x=-2..2);

SymFun(abs(x)^(1/2), x);
plot(abs(x)^(1/2), x=-1..1);
 

@emendes  You can use  makeproc  option. See help on  rsolve  command for this.

@Lali_miani  If you enter in 2d math mode any of  i, j, I  from the Comman Symbols palette, it will work as the imaginary unit. But from the keyboard you have to enter the imaginary unit as  I . I almost never use palettes and only work from the keyboard in 1d math mode (Maple input). 

@minhthien2016  You can make arrows at the ends of the axes of coordinates and labels in the same place if you use  plots:-arrow  and  plots:-textplot  commands.

Example:

Arrow_x:=plots:-arrow([2.3,0],[0.15,0], width=0, head_width=0.12, head_length=0.15):
Arrow_y:=plots:-arrow([0,4.3],[0,0.15], width=0, head_width=0.12, head_length=0.15):
Labels:=plots:-textplot([[2.4,-0.2,"x"],[-0.2,4.4,"y"]], font=[times,16]):
P:=plot(x^2, x=-2..2, color=red, thickness=2, labels=["",""]):
plots:-display(Arrow_x,Arrow_y,P,Labels);

Output:
                         

Edit.