Do a double loop: in the internal one place what you want to suppress, and in the external one everything else and end with a semicolon as in the example:

restart;
for i from 1 to 3 do
for j from 1 to 1 do    
    A := x;
 end;  
B:=C;
end;

It is better to always end a loop with the colon, and if you need to see something in the body of the loop as it is executed, then use  print  as in the example:

for i from 1 to 3 do    
    A := x:
   B := C;
print(%);
end:

 eq:=s=1+2*((tau-t)/T0); 
EQ:=int(f1(t-tau)*(Sum(y[k]*F[k](tau), k = 0 .. M)), tau = t-T0 .. t);
IntegrationTools:-Change(EQ, eq, s);

I do not understand what you are trying to do. Why not just write the following 2 lines and it will be beautiful:

restart;
eq := piecewise(t < 1, sin(t), cos(t));
DocumentTools:-Tabulate(<eq, plot(eq, discont)>, width=50):

Download simpl.mw

Here is a procedure for this:

restart;
MatrixPoly:=proc(nu,M)
local k, Poly, POL; 
for k while k <= M do
Poly[k] := simplify(sum(x^i*GAMMA(nu+1)/(factorial(i)*GAMMA(2*nu)), i = 0 .. k-1));
end do:
POL := <seq(Poly[k], k = 1 .. M)>;
Matrix(M,M, (i,j)->Poly[i]((j-1)/(M-1)));
end proc:

# Examples
MatrixPoly(1,3);
MatrixPoly(1,4);

Edit. Regarding the second problem, it seems I am guessing that you want to create a vector consisting of the values of some function if the argument runs from  a  to  b  through equal steps (total  M  points). If so then

Vec:=(a,b,f,M) -> <seq(f(a+(b-a)*k/(M-1)),k=0..M-1)>:

# Example
Vec(1,2,sin,10)

We can put quotes around function names and omit the square brackets. This will work a little faster:

plot3d('R_ep'(epsilon, t-n, n, 'qN_ep'(epsilon, t-n, n)), n = 1 .. t, epsilon = 10^(-8) .. 10^(-1));

I do not see any bugs in these calculations. Note that  symbolic option works with symbolic expressions. When working with numeric expressions, it is simply ignored. Compare

simplify((p^3)^(1/3),symbolic);
simplify(((-1)^3)^(1/3),symbolic);

                               p
                     1/2+I*sqrt(3)*(1/2)

The second calculation is also correct, as the default Maple works in the complex domain, and returns the principal value of the root (of three values). If you work in the real domain, then use the  surd  function instead of raising to the power 1/3:

surd(-1, 3);

                                 -1

Maple does not simplify the value for (-1)^(1/7)  (in contrast to (-1)^(1/3)), because it is not expressed in radicals, but can be calculated numerically:

simplify((-1)^(1/7));
evalf(%);

                                  (-1)^(1/7)
                    0.9009688678+0.4338837392*I

Maple can also convert   (-1)^(1/7)  to trigonometric form:

convert((-1)^(1/7), trig);
# Or
evalc((-1)^(1/7));

                            cos((1/7)*Pi)+I*sin((1/7)*Pi)
                            cos((1/7)*Pi)+I*sin((1/7)*Pi)

I is protected in Maple (it is the imaginary unit). I replaced I with J :

matrix_inversion_new.mw
 

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