JonesM1 := Compiler:-Compile(
    proc(n::posint)
        local count, num, isPrime, i;

        count := 0;  # Count of primes found
        num := 2;    # Start checking from the first prime number

        while count < n do
            isPrime := true;

            # Check if num is prime
            for i from 2 to floor(sqrt(num)) do
                if num mod i = 0 then
                    isPrime := false;
                    break;
                end if;
            end do;

            # If num is prime, increment count
            if isPrime then
                count := count + 1;
            end if;

            # Move to the next number
            if count < n then
                num := num + 1;
            end if;
        end do;

        return num;  # Return the n-th prime
    end proc
):

JonesM1(1);  # Returns 2
JonesM1(2);  # Returns 3
JonesM1(3);  # Returns 5
JonesM1(4);  # Returns 7
JonesM1(5);  # Returns 11
JonesM1(1000000);

                               2

                               3

                               5

                               7

                               11

                            15485863

n_de_priemgetal_via_compile_en_zonder_28-12-2024.mw

7. The ProgramCode of Bilinear
 with(PDEtools): with(DEtools):
 ## Hirota Bilinear Method
 ## Bilinear Derivative / Hirota Operator
 BD:=proc(FF,DD) local f,g,x,m,opt;
 if nargs=1 then return ‘*‘(FF[]); fi; f,g:=FF[]; x,m:=DD[];
 opt:=args[3..-1]; if m=0 then return procname(FF,opt); fi;
 procname([diff(f,x),g],[x,m-1],opt)-procname([f,diff(g,x)],[x,m-1],opt);
 end:
 ‘print/BD‘:=proc(FF,DD) local f,g,x,m,i; f,g:=FF[];
 f:=cat(f,‘ ‘,g); g:=product(D[args[i][1]]^args[i][2],i=2..nargs);
 if g<>1 then f:=‘‘(g)*‘‘(f); fi; f; end:
 ## collect(expr,f); first!
 getFnumer:=proc(df,f,pow::posint:=1) local i,g,fdenom;
 if type(df,‘+‘) then
 g:=[op(df)];
 fdenom:=map(denom,g);
 for i to nops(fdenom) while fdenom[i]<>f^pow do od;
 if i>nops(fdenom) then lprint(fdenom);
 error "no term(s) or numer=0 when denom=%1",op(0,f)^pow fi;
 g:=numer(g[i]);
 if not type(expand(g),‘+‘) then lprint(g);
 error "Expected more than 1 term about Hirota D-operator" fi;
 return g; fi; lprint(df);
 error "expected 1st argument be type ‘+‘."; end:
 getvarpow:=proc(df::function) local i,f,var,dif,pow; if
 op(0,df)<>diff then lprint(df); error "expected diff function" fi;
 f:=convert(df,D); var:=[op(f)]; dif:=[op(op([0,0],f))];
 pow:=[0$nops(var)]; f:=op(op(0,f))(var[]); for i to nops(var) do
 dif:=selectremove(member,dif,{i});
 pow[i]:=nops(dif[1]);
 dif:=dif[2];
 od; pow:=zip((x,y)->[x,y],var,pow); pow:=remove(has,pow,{0});
 [[f,f],pow[]]; end:
 #convert to Hirota Bilinear Form
 HBF:=proc(df) local i,c,f; if type(df,‘+‘) then
 f:=[op(df)]; return map(procname,f); fi;
 if type(df,‘*‘) then f:=[op(df)];
 f:=selectremove(hasfun,f,diff); c:=f[2]; f:=f[1];
 if nops(f)<>1 then lprint(df); error "need only one diff function factor." fi;
f:=f[]; c:=‘*‘(c[]); f:=getvarpow(f); f:=[c,f]; return f; fi;
 if op(0,df)=diff then f:=getvarpow(df); f:=[1,f]; return f; fi;
 lprint(df); error "unexpected type."; end:
 printHBF:=proc(PL::list) local j,DD,f,C,tmp,gcdC; C:=map2(op,1,PL);
 gcdC:=1; if nops(C)>1 then tmp:=[seq(cat(_Z,i),i=1..nops(C))];
 gcdC:=tmp *~ C; gcdC:=‘+‘(gcdC[]); gcdC:=factor(gcdC);
 tmp:=selectremove(has,gcdC,tmp); gcdC:=tmp[2];
 if gcdC=0 then gcdC:=1 fi; gcdC:=gcdC*content(tmp[1]); fi;
 if gcdC<>1 then C:=C /~ gcdC; fi; DD:=map2(op,2,PL);
 f:=op(0,DD[1][1][1]);
 DD:=map(z->product(D[z[i][1]]^z[i][2],i=2..nops(z)),DD);
 DD:=zip(‘*‘,C,DD); DD:=‘+‘(DD[]); gcdC * ‘‘(DD) * cat(f,‘ ‘,f);
 end:
 ## print Hirota Bilinear Transform
 printHBT:=proc(uf,u,f,i,j,PL,alpha:=1) local DD,g,C,tmp,pl;
 pl:=printHBF(PL); if j>0 then print(u=2*alpha*’diff’(ln(f),x$j));
 else print(u=2*alpha*ln(f)); fi;
 if i>0 then print(’diff’(pl/f^2,x$i)=0); else
 print(pl/f^2=0); fi; NULL; end:
 guessdifforder:=proc(PL::list,x::name)
 local L,minorder,maxorder,tmp; L:=map2(op,2,PL);
 L:=map(z->z[2..-1],L); tmp:=map(z->map2(op,2,z),L);
 tmp:=map(z->‘+‘(z[]),tmp); tmp:=selectremove(type,tmp,even);
 minorder:=0; if nops(tmp[1])<nops(tmp[2]) then minorder:=1 fi;
 tmp:=map(z->select(has,z,{x}),L); tmp:=map(z->map2(op,2,z),tmp); if
 has(tmp,{[]}) then maxorder:=0; else tmp:=map(op,tmp);
 maxorder:=min(tmp[]); fi;
 if type(maxorder-minorder,odd) then maxorder:=maxorder-1 fi;
 [minorder,maxorder]; end:
 guessalpha:=proc(Res,uf,u,f,i,j,PL,alpha) local tmp,res,pl,flag,k;
 flag:=1; tmp:=[op(Res)]; tmp:=map(numer,tmp);
 tmp:=gcd(tmp[1],tmp[-1]); if type(tmp,‘*‘) then
 tmp:=remove(has,tmp,f); fi; if tmp<>0 and has(tmp,{alpha}) then
 tmp:=solve(tmp/alpha^difforder(uf),{alpha});
 if tmp<>NULL and has(tmp,{alpha}) then lprint(tmp);
 for k to nops([tmp]) while flag=1 do
 res:=collect(expand(subs(tmp[k],Res)),f,factor);
 if res=0 then pl:=subs(tmp[k],PL);
 printHBT(uf,u,f,i,j,pl,rhs(tmp[k]));
 flag:=0; fi; od; fi; fi; PL; end:


 Bilinear:=proc(uf,u,f,x,alpha) local su,h,i,j,g1,CB,PL,gdo,DD,Res;
 if hasfun(uf,int) then error "Do not support integral function yet.
Please substitute int function." fi; for j from 0 to 2 do
 Res:=1; su:=u=2*alpha*diff(ln(f),[x$j]);
 h:=collect(expand(dsubs(su,uf)),f,factor);
 if hasfun(h,ln) then next; fi;
 g1:=getFnumer(h,f)/2; g1:=expand(g1); CB:=HBF(g1);
 gdo:=guessdifforder(CB,x);
 for i from gdo[1] by 2 to gdo[2] do
 if i=0 then PL:=CB; else PL:=HBF(int(g1,x$i)); fi;
 DD:=add(PL[i][1]*BD(PL[i][2][]),i=1..nops(PL));
 Res:=collect(expand(diff(DD/f^2,[x$i])-h),f,factor);
 if Res=0 then printHBT(uf,u,f,i,j,PL,alpha); break;
 elif type(alpha,name) and has(DD,alpha) then
 Res:=guessalpha(Res,uf,u,f,i,j,PL,alpha);
 fi; od; if Res=0 then break; fi; od; PL; end:

There is a pdf with examples, so i can test 

wrong code 

Depends on the physical context

restart;
with(Physics):

# 1. Setup and check
# Set default coordinates
Setup(coordinates = [r, theta]);

# Define a metric tensor
Define(g_[~mu, ~nu] = diag(-1, 1, 1, 1)); # Minkowski metric

# Diagnostic test: check a simple case without inert differentials
eq_simple := g_[~mu, ~nu] * v[~mu] * v[~nu]:
SumOverRepeatedIndices(eq_simple); # Expect a contracted result

# 2. Working with inert differentials
# Define the equation with differentials
eq_inert := g_[~mu, ~nu] * Diff(w(r, theta), ~mu) * Diff(w(r, theta), ~nu):

# Diagnostic test: check indices
Check(eq_inert, all);

# Attempt to contract over repeated indices
result_inert := SumOverRepeatedIndices(eq_inert);

# Output the results
print("Result for simple vectors without inertia:");
print(SumOverRepeatedIndices(eq_simple));
print("Result for equation with inert differentials:");
print(result_inert);

# 3. Check various settings
# Attempt explicit metric contraction
Physics:-Simplify(result_inert);

restart;
with(Physics):
Physics:-Setup(inertdiffs = true);
g_[[5, 29, 1]]; 
eq1 := g_[~mu, ~nu] * Diff(w(r, theta), mu) * Diff(w(r, theta), nu);
Physics:-Contract(eq1);
Physics:-Check(eq1, all);

A phaseportrait ...

Maple Tutorial, part 1.2: Phase portrait

with(plots):

# Create the vertical line in segments with labels
L1 := plottools:-line([0, 0], [0, 0.9], color="blue"):
L2 := plottools:-line([0, 1.1], [0, 1.9], color="blue"):
L3 := plottools:-line([0, 2.1], [0, 3], color="blue"):

# Create the disks without transparency
c1 := plottools:-disk([0, 1], 0.1, color="red"):
c2 := plottools:-disk([0, 2], 0.1, color="green"):

# Add labels for each line segment
label1 := plots:-textplot([0.1, 0.45, "Segment 1"], align={left}):
label2 := plots:-textplot([0.1, 1.5, "Segment 2"], align={left}):
label3 := plots:-textplot([0.1, 2.5, "Segment 3"], align={left}):

# Combine the objects and display the plot
display([L1, L2, L3, c1, c2, label1, label2, label3], scaling=constrained, axes=none);

@salim-barzani 
First when you can solve it step by step a problem  , than it should be always possible to make a procedure from it 
Don't know if this makes sense too thi sprocedure 

restart;
with(inttrans):
with(PDEtools):
with(LinearAlgebra):               

undeclare(prime);

with(PDEtools);
declare();                

# Define a procedure that takes the ODE as input
SolveODEWithAdomian := proc(ode)
    local eq, eqs, u0, A, B, j, lap_sol, adomian_sol, final_sol, f, g, m;

    # Step 1: Initial conditions and function definitions
    u0 := (x, t) -> exp(x)*t;  # Define u0 as the base function
    A := Array(0..3);  # Array for Adomian polynomials A[j]
    B := Array(0..3);  # Array for Adomian polynomials B[j]

    # Functions representing the nonlinear terms
    f := u -> u^2;
    g := u -> diff(u, x)^2;

    # Step 2: Calculate the Adomian polynomials A[j] and B[j]
    for j from 0 to 3 do
        A[j] := subs(lambda=0, diff(f(seq(sum(lambda^i * u[i](x, t), i=0..20), m=1..2)), [lambda$j]) / j!);
        B[j] := subs(lambda=0, diff(g(seq(sum(lambda^i * u[i](x, t), i=0..20), m=1..2)), [lambda$j]) / j!);
    end do;

    # Step 3: Perform the Laplace transform of the ODE
    eqs := laplace(ode, t, s);

    # Step 4: Solve in the Laplace domain
    lap_sol := solve({eqs}, {laplace(u(x, t), t, s)});

    # Step 5: Substitute initial conditions for the solution in the Laplace domain
    lap_sol := subs({u(x, 0) = 0, D[2](u)(x, 0) = exp(x)}, lap_sol);

    # Step 6: Transform back to the time domain with the inverse Laplace transform
    adomian_sol := Array(0..3);  # Array to store approximate terms
    for j from 0 to 3 do
        # Define each term of the Adomian solution with inverse Laplace of A[j] and B[j]
        adomian_sol[j] := -invlaplace(1/(s^2) * (A[j] + B[j]), s, t);
    end do;

    # Step 7: Sum all terms of the Adomian solution to obtain the final solution
    final_sol := u0(x, t) + add(adomian_sol[j], j=0..3);

    # Return the result
    return simplify(final_sol);
end proc;

# Define the ODE with the known solution
ode := diff(u(x, t), t $ 2) + u(x, t)^2 - diff(u(x, t), x)^2 = 0;

# Call the SolveODEWithAdomian procedure with the ODE as input
result := SolveODEWithAdomian(ode);

# Print the approximate solution
simplify(result);

wrong example 

# Set parameters
D := 1: # Set an appropriate value for the bending stiffness
rho := 1: # Set an appropriate value for the density
h := 0.1: # Set an appropriate value for the thickness

# Define the differential equation
pde := D*diff(w(x, y, t), x$2) + D*diff(w(x, y, t), y$2) + rho*h*diff(w(x, y, t), t$2) = 0;

# Specify initial and boundary conditions
ic := w(x, y, 0) = sin(Pi*x)*sin(Pi*y), D[1, t](w(x, y, t)) = 0; # initial displacement and velocity
bc := w(0, y, t) = 0, w(1, y, t) = 0, w(x, 0, t) = 0, w(x, 1, t) = 0; # boundary conditions

# Solve and visualize
sol := pdsolve([pde, ic, bc], w(x, y, t), numeric);

# Plot the vibration
plots:-animate3d(subs(sol, w(x, y, t)), x = 0 .. 1, y = 0 .. 1, t = 0 .. 1, frames = 30);

As example , you can get all sorts of functions.

f := z -> z * hypergeom([1, 1], [2], -z);

f := proc (z) options operator, arrow; z*hypergeom([1, 1], [2], 

   -z) end proc


f := z -> z * hypergeom([1, 1], [2], -z):
simplify(f(z));

                           ln(1 + z)

Now as procedure

restart;

solveODE := proc(data)
    local the_result, ode, current_ode, sol, item;
    
    # Create an empty array to store results
    the_result := Array(1..0):
    # Header for the table
    the_result ,= [P, Q, R, "current ode", "solution"]:
    
    # Differential equation
    ode := diff(F(xi), xi)^2 = P * F(xi)^4 + Q * F(xi)^2 + R:
    
    # Iterate through each set of P, Q, R in data
    for item in data do
        current_ode := limit(eval(ode, item), m = 1, 'left'):
        sol := [dsolve(current_ode)]:
        the_result ,= [rhs(item[1]), rhs(item[2]), rhs(item[3]), current_ode, sol]:        
    od:
    
    # Convert the results to a list and display as a table
    the_result := convert(the_result, listlist):
    return DocumentTools:-Tabulate(the_result, weights = [10, 20, 20, 45, 45], width = 60):
end proc:

# Run the procedure with the given data
data := [
    [P = m^2, Q = -(1 + m^2), R = 1],
    [P = 1, Q = -(1 + m^2), R = m^2],
    [P = -m^2, Q = 2 * m^2 - 1, R = 1 - m^2],
    [P = 1, Q = 2 * m^2 - 1, R = -m^2 * (1 - m^2)],
    [P = 1/4, Q = (m^2 - 2) / 2, R = m^2 / 4],
    [P = m^2 / 4, Q = (m^2 - 2) / 2, R = m^2 / 4]
]:

solveODE(data);

I have here an example of a table with a complex function, in which the Explore command is inserted in the table cell.
An absolute complex plot is not attached 
It's more of a draft worksheet to get some ideas 

mprimes_parameter_13-10-2024-salim.mw

@Susana30 
Yes, you can copy it from here. 
To get floating point numbers( decimal)  , add a  period with zero ( i did for  2 ) for one intersection point  in Cartesian coordinates.
 

theta_values := [2.0*Pi/3, 4*Pi/3]:  # The two theta values for the intersection points

with(plots):

# Exercise: Calculate the intersection points in Cartesian coordinates
# Step 1: Solve the equation to find theta
theta_values := [2.0*Pi/3, 4*Pi/3]:  # The two theta values for the intersection points

# Step 2: Calculate the radius values (r) for the found angles
r_circle := map(t -> -6*cos(t), theta_values):  # For the circle

# Step 3: Calculate the Cartesian coordinates for the two intersection points
points_cartesian := [
    [r_circle[1]*cos(theta_values[1]), r_circle[1]*sin(theta_values[1])],  # Intersection point 1
    [r_circle[2]*cos(theta_values[2]), r_circle[2]*sin(theta_values[2])]   # Intersection point 2
]:

# Display the calculated intersection points in Cartesian coordinates
points_cartesian;  # This will display the two points as (x, y)

# Expected output:
# [-1.5, 2.598], [-1.5, -2.598]

# Step 4: Plot the circle and cardioid in polar coordinates
p2 := polarplot([-6*cos(theta), 2 - 2*cos(theta)], theta = 0..2*Pi, color = [green, purple], legend = ["Circle", "Cardioid"]):

# Step 5: Plot the intersection points in Cartesian coordinates
intersection_points := pointplot(points_cartesian, symbol=circle, color=red, symbolsize=15):

# Step 6: Combine the plots
display(p2, intersection_points);

Download matrix_sommen_-maple_primes.mw