This can be done in many ways. Here is a solution using a procedure that I wrote many years ago. All comments in the procedure are written in Russian. The full text of the procedure itself can be found in the attached file, and here its application to your example:

Download QuadricCurveAnalysis.mw

This is a known bug in  RealDomain  package when it does not take into account the domain of the equation. Here is a simple workaround:

restart;
RealDomain:-solve({(x-1)*sqrt(x^2-4) = 0, x^2-4>=0}, x);

                                     {x = -2}, {x = 2}

restart;
f:=t->piecewise(t<1 and t>=-1,t/2+1/2,t>=1 and t<2,-t+2);
F:=t->f(t-floor((t+1)/3)*3); #  F  is a periodic extension of  f  to the whole real axis
F(4.5);
plot(F(t), t=-4..5, scaling=constrained, size=[1200,150]);

Use the  LinearAlgebra:-GenerateMatrix  command for this (3 steps):
 

Download MatrixForm.mw

If you do not like the order of the elements in the set, then you can convert it into the list, and then sort this list according to the desired order:

E:={{1,2,3,4,5,6},{4,5,6,7,8,9,10,11},{10,11,12,13,14}};
        E := {{10, 11, 12, 13, 14}, {1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}}
E1:=convert(E, list);
       E1 := [{10, 11, 12, 13, 14}, {1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}]
sort(E1, (x,y)->x[1]<y[1]);
       [{1, 2, 3, 4, 5, 6}, {4, 5, 6, 7, 8, 9, 10, 11}, {10, 11, 12, 13, 14}]

Of course, to avoid this hassle, you should immediately specify E as a list and not as a set.

See help on the  diff  command.

The code for your example:

w:=exp(-beta*t-mu*x)*u(x,t);
diff(w, t);

Eq:=v1(t) - R1*(-i3(t)/n13 - i2(t)/n12) - L1*diff(-i3(t)/n13 - i2(t)/n12, t) = n12*(v2(t) - R2*i2(t) - L2*diff(i2(t), t));
Eq1:=expand((lhs-rhs)(Eq));
Terms:=[op(Eq1)];
DTerms:=select(has, Terms, diff);
Lhs:=add(DTerms);
collect(Lhs, diff)=collect(-(Eq1-Lhs), R1);  # The final result

You must specify the integration range, but you cannot use the symbol to as a name because it is a reserved word. I replaced it with  a :

a:=3:
evalf(int(exp((5+I*6)*x)*sin(1+erf(x)),x=-a..a));
evalf(int(exp((5+I*6)*x)*cos(1+erf(x)),x=-a..a));

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
subs(e2=e1/n12, A);

Or more automatically:

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
eqn2 := n12 = e1/e2;
subs(e2=solve(eqn2, e2), A);

Or using the  applyop  command:

restart;
eqn1 := e1*i1 + e2*i2 + e3*i3 = 0;
A:=i1=expand(solve(eqn1, i1));
applyop(simplify, [2,1], A, {e1/e2=n12});

Edit.
 

Use  plots:-pointplot3d  and  plots:-display  commands for this.

An example:

pointplot3d.mw

You forgot to specify values for  xbas, tbas :

Download anim_plot3d.mw
 

I rewrote your code as 2 procedures. The first one generates a full random polynomial of degree  N  with coefficients from the range  R . The second finds all rational roots of a polynomial with integer coefficients:

Download rational_zeros_new.mw

We can easily solve this problem using the  seq  command with the  `if` - condition:

L:=[$ 1..50];
L1:=[seq(`if`(irem(i,2)=0,NULL, L[i]), i=1..nops(L))];
L2:=[seq(`if`(irem(i,3)=0,NULL, L1[i]), i=1..nops(L1))];

Edit. This method is about 2 times faster than using the  remove  command:

L:=[$ 1..1000000]:
CodeTools:-Usage([seq(`if`(irem(i,2)=0,NULL, L[i]), i=1..nops(L))]):
CodeTools:-Usage(RemoveEveryNth(L, 2)):

I have simplified your procedures somewhat by removing unused formal parameters and unnecessary local variables. Note also that the external call of  packages does not work in the procedures. Use  the  uses  command to call the package inside the procedure.

Download func_new.mw

1. First, you made a common syntax error: the exponent is encoded in Maple as  exp(x)  not  e^x  (in Maple  e  is only a symbol rather than the number  2.71828...).

2. Your integral contains 4 parameters (a, b, c, d), but their role is different. The most important is the parameter  b and we can easily get a solution (in many cases) depending on the value of  this parameter. The examples below show that the solution form depends heavily on the value of  b :

restart;
f:=b->d*exp(-b*x)/(((a*exp(-2*b*x)+c*exp(-4*x))));
F:=b->int(f(b),x);

 # Examples:
for b from 0 to 5 do
F(b);
od; 

Edit.