Several points to help you along:

1. The output of Pollard's p-1 algorithm, when it works, is simply a proper factor of n. It may be composite or prime.
2. The algorithm only works if n has a prime factor p such that p-1 has only relatively small prime factors. No practical cryptography system uses such n. "Relatively small" means less than B as used in my procedure below. If B is set to n^(1/6), the probability that the algorithm will work is 1/27 (see the Wikipedia article).
3. Detecting primality is much easier than finding a factor. The algorithm shouldn't be used to test primality. The algorithm should reject input that is prime.
4. Since you use 2 as the base of your exponentations, your procedure needs to reject even input.
5. To find the gcd of integers, use igcd rather than gcd. The latter will work, but it is slower as it is intended for polynomials.
6. To do a modular exponentiation with a large exponent, use a &^ b mod n rather than a^b mod n. The former is much faster.

Here's my Pollard procedure, based on the Wikipedia article Pollard's p-1 algorithm:

Pollard:= proc(
     n::And(posint, odd, Not({identical(1), prime})),
     {B::posint:= max(trunc(n^(1/4)), 10^4)}
)
description "Pollard's p-1 algorithm";
local q:= 1, lnn:= evalf(ln(n)), g, a:= 2;
     while q <= B do
          q:= nextprime(q);
          a:= a &^ (q^trunc(lnn/ln(q))) mod n;
          g:= igcd(a-1,n);
          if g > 1 then return g end if
     end do;
     FAIL
end proc;

Do you mean this?

plots:-display(
     plots:-intersectplot(x+y=1, x^2+y^2+z^2=1, x= -1..1, y= -1..1, z= -1..1),
     plot3d(1, theta= 0..2*Pi, phi= 0..Pi, coords= spherical),
     scaling= constrained, axes= boxed
);

Simply use coeff. There's nothing special about it being a Lie derivative; it's treated as any other sum of terms.

coeff(L1fh, D_x);

You should use explicit multiplication signs and not rely on juxtaposition implying multiplication.

int(1/(x*(1-x))*(1+2*x), x);

1/2*u^2*(1-u)^2*(1+2*u)^2;

Now the integration is performed, and the second expression is interpretted correctly.

teksbiasa:= `Hello! Bob`:
ListTools:-Flatten(
     ListTools:-Pad[1](
          (w-> w+~nops(w))~(convert~(StringTools:-Split(teksbiasa, ` `), bytes)),
          32
     )[..-2]
);

The following plots the first 5122 zeros (the plot in the paper shows 10,000) in about 18 minutes. To get more zeros, expand the imaginary range, which I set at -9999..9999. The zeros are rather evenly distributed throughout the imaginary range.

R:= [RootFinding:-Analytic(add(1/n^s, n= 1..5), s, re= -3..1, im= -9999..9999)]:
nops(R);

     5122

To normalize the roots modulo 2*Pi*I/ln(5), as is done in the paper, use frem.

M:= evalf(2*Pi/ln(5)):  
Frem:= proc(x,y) local r:= frem(x,y); `if`(r<0, r+y, r) end proc:

One way to plot complex points (the easiest way as far as I'm concerned) is to separate them into real and imaginary parts with Re and Im.

P:= map(z-> [Re(z), Frem(Im(z),M)], R):
plot(P, style= point, symbol= point, symbolsize= 1);

Sure, it's easy:

sumif:= (L::{list,set}, B)-> `+`(select(B, L, _rest)[]):

Now let's say that I have a list:

R:= rand(9):  L:= ['R()' $ 9];

And I want to sum all those elements that are greater than 5:

sumif(L, `>`, 5);

 30

Another way: If L is the list of digits, do

(parse@cat@op)(L);

ListTools:-Reverse(convert(mylist, base, 10));

The reverse of [2,7,5] is [5,7,2]. You have it as [5,7,5].

Why is 3 any different than the 11 that I showed you how to add to every term last week?

nilaiASCII +~ kekuncirahsia;

What you want can be done with a single simplify with side relations command. First you need to freeze the expressions e^epsilon and e^(-epsilon) so that every expression in the matrix is polynomial. (Note that these aren't the same as exp(epsilon) and exp(-epsilon); however, I've continued with your erroneous usage.) The expression that you have labeled (3) I call ex. Then all that you need to do is

M1:= subs(e^epsilon= freeze(e^epsilon), e^(-epsilon)= freeze(e^(-epsilon)), M):
thaw(simplify(ex, {seq(seq(T[i]*U[j]= M1[i,j], i= 1..op([1,1], M)), j= 1..op([1,2], M))}));

Since theta(t) only appears in derivatives, it is easy to reduce the order. Change diff(theta(t), t) to dtheta(t) and diff(theta(t), t, t) to diff(dtheta(t), t):

eq1:=J*diff(dtheta(t),t)+b*dtheta(t)=K*i(t):
eq2:=L*diff(i(t),t)+R*i(t)=V(t)-K*dtheta(t):
DCMotor:=[eq1,eq2];
Sys:=DiffEquation(DCMotor,[V(t)],[dtheta(t)]);

The following command will take care of the asymptotes automatically:

Student:-Precalculus:-RationalFunctionPlot((x-1)/(x-2));

alpha__0:= x-> sqrt(d/x)/2*exp((-D__s - i)*Omega*x/V__La):
alpha__90:= x-> 3*sqrt(d/32/x)*exp((-D__s - i)*Omega*x/V__s):
alpha__theta:= (x,theta)-> alpha__0(x)*cos(theta)^2 + alpha__90(x)*sin(theta)^2:

If you intend for the i to refer to the complex unit sqrt(-1), then you'll need to change it to I.