Here's your formula in Maple. Note that the angle phi is measured in degrees.

Fi:= proc(phi)
local p:= phi*Pi/180, b:= 0.081819221, sb:= b*sin(p);
     7915.7*log10(tan((Pi/4+p/2)*((1-sb)/(1+sb))^(b/2)))
end proc;

I think that the numeric integration is getting hung up trying to achieve high accuracy, which you don't need for a plot. I changed the last two lines of your code to

dpp:= Int(Re(dpdx), x= 0..1, digits= 4, method= _d01ajc):
dpp2:= eval(dpp, K= -0.1);
plot(dpp2, Q= -1..1, axes= box,color= [blue]);

and I got the plot in under ten seconds. Not necessarily all of the changes that I made are necessary for getting a fast plot.

The evalm command has been obsolete for over a decade and has been deprecated for years. If you show the code that you are trying to use, I can update it.

Here is my implementation of Rodrigues's formula. K is the axis of rotation.

restart:
Rodrigues:= proc(K::Vector(3),v::Vector(3),theta::algebraic)
local
     k:= LinearAlgebra:-Normalize(K,2),
     kx:= < < 0, k[3], -k[2] > | < -k[3], 0, k[1] > | < k[2], -k[1], 0> >
;
     v + kx.v*sin(theta) + (kx^2).v*(1-cos(theta))
end proc:
           
R:= unapply(Rodrigues(<1,1,1>, <1,1,0>, t), t):
plots:-spacecurve(R(t), t= -Pi..Pi);

Your system is autonomous, and it is only differentiated with respect to the first variable. It is thus equivalent to a system of ODEs. The ODEs can be easily solved.

sys:= {eq||(1..3)}:
sysA:= subs([a*t= NULL, a^2*t= NULL], evalf(sys)):
dsolve(sysA, convert_to_exact= false);

I don't know how to procede from this point. It is not even clear to me that there is any way to procede.

Use convert(%, StandardFunctions) followed by simplify(%) to obtain

which is equivalent to the desired form that you posted.

Here's a procedure to generate it for an arbitrary Vector or list of xs.

B:= proc(X::{Vector,list})
local n:= numelems(X);
     LinearAlgebra:-Determinant(
          Matrix(
               n,
               (i,j)-> `if`(j < i, -1, binomial(n-i,j-i)*X[j-i+1]),
               shape= Hessenberg[upper]
          )
      )
end proc:

B(Vector(4, symbol= x));

dsolve({diff(y(x),x)=x*y(x)^2, y(1)=2});

taylor(rhs(%), x=1);

Here's an idea to deal with this with a modicum of generality:

restart:
Int(1/x^n, x= 10..100);

A:= value(%):
piecewise(seq([n = d, limit(A, n= d)][], d= discont(A,n)), A);

Are you just interested in real roots? Looking at plot([tan(x/2), tanh(x/2)]) it is clear that there is exactly one real root for every period of tan(x/2). So there is exactly one real root in each interval (2*k-1)*Pi..(2*k+1)*Pi for any integer k. Once you decide on the interval, the root can be found with fsolve:

fsolve(eq, x= 5*Pi..7*Pi);

20.42035225

If you want the complex roots, let me know. It's not as easy to analyze.

I added the plot of the arc, c2, to your two plots, and specified that it be filled.

c1:= plottools:-circle([0, 0], 1, color = red):
c2:= plot([cos(t), sin(t), t= -arccos(1/2)..arccos(1/2)], filled):
p2:= plots:-implicitplot(x = 1/2, x = -2 .. 2, y = -1 .. 1.1, colour = blue, linestyle = 1, thickness = 2);
plots:-display([c1, c2, p2]);

The following procedure will change the _Cs for you, regardless of how many there are.

ChangeC:= proc(Sol, {oldprefix::symbol:= _C, newprefix::symbol:= C})
local n:= length(oldprefix);
     subsindets(
          Sol,
          And(
              symbol,
              satisfies(
                   proc(x::symbol)
                       try
                            substring(x, 1..n) = oldprefix
                            and parse(substring(x, n+1..-1))::nonnegint
                       catch:
                            false
                       end try
                   end proc
              )
          ),
          x-> newprefix[parse(substring(x, n+1..-1))]
     )
end proc:

Examples:

Sol:= dsolve(diff(y(x),x$2)=y(x)):
ChangeC(Sol);

ChangeC(Sol, newprefix= A);

What I would do is place the less relevant sections of code in Code Edit Regions. These can be shrunk to arbitrarily small size while still remaining scrollable. You can find Code Edit Region on the Insert menu. Remember to end statements in the Regions with a colon to supress their output. The output, if any, appears after the Region.

Remove the storage= sparse from your Matrix declaration. Apparently sparse matrices do not work with plots:-pointplot. Anyway, your Matrix is completely filled in, so there's no point in using sparse storage.

By "as a series" do you mean side-by-side plots? That's requires a slight change to Mehdi's Answer:

plots:-display(< F1 | F2 | F3 >);

I can't upload these array plots.

Listing all 3.6 million elements of S10 is silly. Listing (or even storing in memory) all ~10^18 elements of S20 is impossible.

If you want to list all permutations on n elements (for reasonably small n), you can do it by

S||n:= combinat:-permute(n);

Each of these permuations can be converted to disjoint cycle notation by convert(..., disjcyc):

map(convert, S||n, disjcyc);

A more reasonable goal is to create structures in Maple that Maple understands to be the groups S10 or S20 without creating every member. You can do this with permgroup. The symmetric group on n elements (S||n) can be generated by two permutations: a rotation (cycle) of all the elements and a transposition of any two.

S10:= permgroup(10, {rot= [[$1..10]], trans= [[1,2]]});

This group object and its members can be manipulated by the commands in package group. See ?group . (Oddly, permgroup is a top-level command, not explicitly part of the package.)

Could you be more specific about making a list of transpositions and cycles?