You can use the Student:-Calculus1:-Roots command for this:

restart;
f(x) := 2.0*cos(1.5*x + 4.0) + 2.2:
g(x) := 2.0*cos(1.6*x + 3.85) + 2.0:
Student:-Calculus1:-Roots(f(x)=g(x), x=-1..4);             

                                                    [-0.7967442905, 3.834531623]

In Maple we can use the integral test (see wiki  https://en.wikipedia.org/wiki/Convergence_tests ):

int(1/k^(2-cos(1/k)), k=1..infinity);

In the list of your equations, you missed  eqn8 . Additionally, you can use the  fsolve  command rather than the  solve  one :

restart;
eqn1 := W__1 + W__2 + W__3 + W__4 = 4;            
eqn2 := W__1*zeta__1 + W__2*zeta__2 + W__3*zeta__3 + W__4*zeta__4 = 0;
eqn3 := W__1*zeta__1^2 + W__2*zeta__2^2 + W__3*zeta__3^2 + W__4*zeta__4^2 = 2/3;
eqn4 := W__1*zeta__1^3 + W__2*zeta__2^3 + W__3*zeta__3^3 + W__4*zeta__4^3 = 0;
eqn5 := W__1*zeta__1^4 + W__2*zeta__2^4 + W__3*zeta__3^4 + W__4*zeta__4^4 = 2/5;
eqn6 := W__1*zeta__1^5 + W__2*zeta__2^5 + W__3*zeta__3^5 + W__4*zeta__4^5 = 0;
eqn7 := W__1*zeta__1^6 + W__2*zeta__2^6 + W__3*zeta__3^6 + W__4*zeta__4^6 = 2/7;
eqn8 := W__1*zeta__1^7 + W__2*zeta__2^7 + W__3*zeta__3^7 + W__4*zeta__4^7 = 0;

fsolve({eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8}, {W__1, W__2, W__3, W__4, zeta__1, zeta__2, zeta__3, zeta__4});

  {W__1 = 1.620164810, W__2 = 0.3798351902, W__3 = 1.620164810, W__4 = 0.3798351902, zeta__1 = -0.1911644819, zeta__2 = 0.8495280449, zeta__3 = 0.1911644819, zeta__4 = -0.8495280449}

Edit.  You can also get exact symbolic solutions by using the  explicit option. But there will be 24 such solutions:

Sys:={eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8}:
solve(Sys, explicit);
evalf([%]);
nops(%);

Always use the add command when summing a finite number of specific terms. The  sum  command is usually used to sum an infinite number of terms or a finite, but not predetermined (a symbolic summation).

restart;
U[0] := x^2;
for k from 0 to 1 do U[k+1] := add(U[s]*(diff(U[k-s], x)), s = 0 .. k) end do;

# U[0]*(diff(U[0], x));
                            
# U[1];

Addition. Below are 2 simple examples where the sum  command is needed:

restart;
sum(1/k^4, k=1..infinity);
sum(k^4, k=1..n);
simplify(%);

It seems the method _d01ajc is only suitable for real-valued functions. Remove this option and the error disappears. In your example, some logarithms will be from negative numbers. Below is a simple example with a similar error:

evalf(Int(ln(x-2), x=0..1, method=_d01ajc));
evalf(Int(ln(x-2), x=0..1));

Error, (in evalf/int) unable to obtain a real result
                  0.3862943611 + 3.141592654 I
 

We have to give names for each root. Replace the last line of your code with the following line:

assign(seq(z[i] = sol1[i], i = 1 .. 4));

You can do this programmatically as in the example below:

restart;
Sys:={x+y-z = 3, x-y-z = 5, -x-y-z = 7};
%piecewise(``, Sys[1], ``, Sys[2], ``, Sys[3]);

In the example below, we first generate a list of 30 random numbers, each in the range from 0 to 5, and then select a sublist of numbers from 0 to 1 from it, indicating the indices:
 

restart;
r:=rand(0...5.):
y:=[seq(r(), n=1..30)];
[seq(`if`(y[n]>0 and y[n]<1,[n,y[n]],NULL), n=1..nops(y))];

  y := [0.2499109628, 2.078923730, 2.252410250, 0.3475034102, 

    1.766288344, 3.629756492, 3.635102067, 3.892007975, 

    3.320931315, 0.008823971507, 3.465092938, 1.209628080, 

    1.722913437, 3.215511858, 4.467963158, 0.5108327385, 

    0.2750667362, 3.363685982, 2.075595849, 1.224136121, 

    3.783826445, 0.5091763968, 3.694547905, 0.2999657008, 

    2.685906488, 1.617331965, 0.2915647810, 1.493696599, 

    4.158416643, 3.518466880]
 [[1, 0.2499109628], [4, 0.3475034102], [10, 0.008823971507], 

   [16, 0.5108327385], [17, 0.2750667362], [22, 0.5091763968], 

   [24, 0.2999657008], [27, 0.2915647810]]
 

restart;
with(Statistics):
x := RandomVariable(Binomial(45, 0.9)):
for i from 0 to 45 do
    P[i]:=[i,ProbabilityFunction(x, i)];
end do:
Points:=convert(P, list);
plot(Points, style=pointline, size=[800,400], labels=["i","x"],labelfont=[times,16]);

   
Edit.

You don't need the combinat package for this. You can just use a nested loop for this (triple loop for your example since k = 3):

restart;
L:=[1,2,3,4,a]: n:=nops(L):
for i1 from 1 to n-2 do
for i2 from i1+1 to n-1 do
for i3 from i2+1 to n do
print([L[i1],L[i2],L[i3]]);
od: od: od:

                           [1, 2, 3]
                           [1, 2, 4]
                           [1, 2, a]
                           [1, 3, 4]
                           [1, 3, a]
                           [1, 4, a]
                           [2, 3, 4]
                           [2, 3, a]
                           [2, 4, a]
                           [3, 4, a]
 

Try the following. In Maple 2018, both methods lead to the same result:

restart;
f := x*y:
JH := int(Heaviside(f-1/2), [x=0..1, y=0..1]); 
JP := int(piecewise(x*y>1/2, 1, 0), [x=0..1, y=0..1]);

restart;
y1:=2*x-4: y2:=-2*x+4:
A:=plot([y1,y2], x=0..4, -2..2, color=blue, thickness=2, labels=[``,``]):
B:=plottools:-sector([2, 0], 0.5, arctan(2) .. Pi-arctan(2), color = pink):
T:=plots:-textplot([[4.3,-0.2,x],[-0.2,2.3,y]], font=[times,16]):
plots:-display(A,B,T, view=[-0.2..4.3,-2..2.3], scaling=constrained);

Or

min(op(xLp_tmp));

restart;
simplify(-0.27059805007310 * sin(0.12 + epsilon) +0.27059805007310 * sqrt(1.-cos(0.12 + epsilon) ^ 2)) assuming sin(0.12 + epsilon)>=0;

                                                           0.

or

expr:=-0.27059805007310 * sin(0.12 + epsilon) +0.27059805007310 * sqrt(1.-cos(0.12 + epsilon)^2);
subs(0.12 + epsilon=t, expr);
simplify(%) assuming t>=0 and t<=Pi;

restart;
G:=GraphTheory:-RandomGraphs:-RandomGraph([seq(v[i],i=1..8)], 10, connected);
GraphTheory:-DrawGraph(G);

Edited.