Use  factor  instead of  collect :

A:=factor(c4*dnub*kpbr*ksr*nur*nurdel + c4*dnur*kpbr*ksr*nub*nurdel);
`*`(op(1..3,A))*expand(`*`(op(4..5,A)));

                       A := c4*kpbr*ksr*nurdel*(dnub*nur+dnur*nub)
                      c4*kpbr*ksr*(dnub*nur*nurdel+dnur*nub*nurdel)

restart;
f:=(a+b)*x1+(a^2+b^2)*(x1+x2^2)+c*(x2-x3)*a*b:
p:=[a+b,a^2+b^2,c*b*a]:
R:=p=~1:
cof:=[seq(algsubs(R[i], [op(f)][i]),i=1..nops(p))];

                    cof := [x1, x2^2+x1, x2-x3]

In Maple 2018 we can make it significantly shorter (I mean vv's answer):

Op := cosh(cos(2*t)); 
simplify(convert(Op, exp)); 
map(int, %, t = 0 .. 2*Pi);

Edit. The same code works in Maple 2015 - 2017. It probably also works in Maple 2019 (I don't have this version).

You have to give different names for different values of  k , so I replaced  Sol_f[H]  with  Sol_f[k] . Next, you should get the data for each  k  in the form of a Matrix or Array. In the example, when creating the matrix  M , I divided the range  x=0..1  into 10 parts. Of course you yourself can choose any number of parts. Below I show how this can be done for k = 1. You can export these data to Excel using  ExcelTools:-Export  command. Read help for this.

Download Untitled_new.mw

Try:
<x^`1`, x^`2`, x^`3`>;

or it can be done programmatically if the vector is already created:

X:=<x^1,x^2,x^3>;
map(t->`if`(degree(t)=1,t^`1`,op(1,t)^convert(op(2,t), symbol)), X);

Edit.

Use the  iterationlimit  option.
An example  (the first 10 iteration):

We see that already at the 8th iteration, the desired accuracy was achieved.

Edit.

Download number_iter.mw

with(plots):
with(plottools):
setoptions(style=patch,view=[-5..5,-5..5], scaling=constrained):
basis_i:=polygon([[0,0],[1,0],[0,1]],color=red):

N:=20:
F:=(x,y) -> (1-k/N)*<x,y> + k/N*(A.<x,y>):
L := transform((x,y)->convert(F(x,y),list)):
A:=<<2,1>|<-1,1>>;
                          
frames:=seq(display(L(basis_i)),k=0..N):

display(frames, insequence=true);

Or

with(plots):
with(plottools):
setoptions(style=patch,view=[-5..5,-5..5], scaling=constrained):
basis_i:=polygon([[0,0],[1,0],[1,1],[0,1]],color=red):

N:=20:
F:=(x,y) -> (1-k/N)*<x,y> + k/N*(A.<x,y>):
L := transform((x,y)->convert(F(x,y),list)):
A:=<<2,1>|<-1,1>>;
                          
frames:=seq(display(L(basis_i)),k=0..N):

display(frames, insequence=true);

Or here is another technique for animating linear transformation of a square grid:

restart;
with(plots):
with(plottools):
Square_grid:=display(seq(line([0,i],[5,i],color=red),i=0..5),seq(line([j,0],[j,5],color=red),j=0..5)):

F:=(x,y) -> (1-t)*<x,y> + t*(A.<x,y>):
L := transform((x,y)->convert(F(x,y),list)):
A:=<<2,1>|<-1,1>>:

animate(display,[L(Square_grid)], t=0..1, frames=60, scaling=constrained, size=[700,500]);

If you want to get an animation of a n-th iteration of this mapping, replace the last line of code with

animate(display,[(L@@n)(Square_grid)], t=0..1, frames=60, scaling=constrained, size=[700,500]);

where  n>=2  is an integer.

Edit.

Edit.

Download stat1_new1.mw

For a better understanding of the rule of the numbering of these points, I made an animation of this process. I reduced the number of points to 70, otherwise the file turned out to be too large and there were problems. We see that the numbering occurs in a spiral of Archimedes, which spins counterclockwise with an angular pitch of approximately 138 degrees:

sunflower_new1.mw

 Edit.    

We use polar coordinates. First we find  r0  corresponding to the intersection of these surfaces:

Download polar.mw

The logic is very simple. The   plots:-implicitplot  command works by calculating the values of a function of two variables  F(x,y)  on a grid. For example, if  F(x1,y1)<0  and  F(x2,y2)>0 , then on the segment connecting these points there exists a point such that  F(x0,y0)=0  (of course we assume that the function  F  is continuous.
Your function  F(r,theta)=r-cos(theta)  in polar coordinates . It is equal to 0 at points of the circle with center  (x,y)=(1/2, 0)  and radius 1/2. You are plotting an region outside this circle. I replaced  r=1/3  in your code with  r=1.5  so that the whole circle is visible. For unknown reason, Maple does not paint over the region to the left of the Oy axis. For the correct plotting, we can use the Cartesian coordinates:

Addition. If there is no sign change, then the  plots:-implicitplot  command fails. Compare (obviously abs(x-y)=0 is equivalent to x-y=0) .

plots:-implicitplot(abs(x-y)=0, x=-1..1,y=-1..1);
plots:-implicitplot(x-y=0, x=-1..1,y=-1..1);

Download region_polar.mw

An easy way is to write both equations in one system. Here is an example:

Edit. If you already have separate solutions for these equations (the same equations), then you can do so

restart; 	
Sol1:=dsolve({diff(P(x),x)=2*x,P(0)=0},numeric):
Sol2:=dsolve({diff(V(x),x)=cos(x),V(0)=0},numeric):
plot([t->eval(P(x),Sol1(t)), t->eval(V(x),Sol2(t)),0..10]);

The result is the same as above.

Download PV.mw

We can make the animation much smoother if we provide linear interpolation between successive frames:

restart;     
SOLNSuy[1, 1] := 2.5872902469406659197*10^(-20)-.65694549571241255901*y
+1.9708364871372376767*y^2-1.3138909914248251176*y^3-1.6010739356637904911*10^(-19)*y^4:     
SOLNSuy[2, 1] := -4.002204462000*10^(-20)-1.7879176897079605225*y
+5.3637530691192141414*y^2-3.5758353794044226250*y^3-6.8309939211286845440*10^(-12)*y^4:     
SOLNSuy[3, 1] := -1.1953264450000*10^(-19)-3.2481690589079594122*y
+9.7445071767154794599*y^2-6.4963381177952273213*y^3-1.2292726248071398400*10^(-11)*y^4:    
SOLNSuy[4, 1] := -2.6720465500000*10^(-19)-4.9239979672954025921*y
+14.771993901873204315*y^2-9.8479959345587718955*y^3-1.9029826928878336000*10^(-11)*y^4:     
SOLNSuy[5, 1] := 3.416928541000*10^(-20)-6.7268498492441931137*y
+20.180549547714413714*y^2-13.453699698443639810*y^3-2.6580790570532587008*10^(-11)*y^4:     
SOLNSuy[6, 1] := -2.554122292000*10^(-20)-8.5884528335125514887*y
+25.765358500514014457*y^2-17.176905666966875698*y^3-3.4587270427710613504*10^(-11)*y^4:     
SOLNSuy[7, 1] := -9.206107680000*10^(-20)-10.456823708331499352*y
+31.370471124965259849*y^2-20.913647416590986491*y^3-4.2774005353527132160*10^(-11)*y^4:     
SOLNSuy[8, 1] := 1.9644186790000*10^(-19)-12.293003938471349390*y
+36.879011815379230436*y^2-24.586007876856948223*y^3-5.0932823222176363520*10^(-11)*y^4:    
SOLNSuy[9, 1] := -3.775112769000*10^(-19)-14.068404975282556550*y
+42.205214925807397100*y^2-28.136809950465931724*y^3-5.8908824448577377280*10^(-11)*y^4:     
SOLNSuy[10, 1] := 1.146281780000*10^(-19)-15.762658869974768890*y
+47.287976609878780960*y^2-31.525317739837422477*y^3-6.6589592851037286400*10^(-11)*y^4:

t:=a-floor(a):
F:=a->(1-t)*SOLNSuy[floor(a), 1]+t*SOLNSuy[floor(a)+1, 1]:
     
plots[animate](plot, [F(a), y=0..1], a = 1 .. 10, frames=50, size=[800,400]);

If it is necessary to introduce a new function (let's say with the name  g )  such that  g(zeta)=f(x)  under the condition  zeta=x/a , then it is obvious that  g:=zeta->f(a*zeta)
So the new function  g  will be

f:=x->A1*sin(k*x)+A2*cos(k*x)+A3*sinh(k*x)+A4*cosh(k*x):
g:=unapply(f(a*zeta), zeta);

If we need to select terms of a certain degree (relative to some variable) in a polynomial, then this can be done as in this example:

P:= y^2 - x^2/y:
n:=degree(P,x);
select(t->degree(t,x)=n, P);