I'm guessing that the original expression is missing the parentheses (for the numerator and denominator), as well as the multiplication signs. Using the indication in the problem statement, we first find the tangency points, and then the tangent lines at these points:

restart;
f:=x-> (x^3+9*x^2-9*x-1)/(x^4+1):
solve(D(f)(x)=0.5);
L:=select(type,[%], realcons); # the tangency points
Lines:=map(p->f(p)+0.5*(x-p), L); # the tangency lines

plot([f(x),Lines[]], x=-6..6,-4..10, color=[red,blue,green,brown,cyan], scaling=constrained, legend=[f(x),Lines[]]); 

0.4396034712, 1.545359967, 1.135924545 + 3.219224078 I, 0.07028667536 + 0.8729569266 I, -0.8928416048, -3.504544275,  0.07028667536 - 0.8729569266 I, 1.135924545 - 3.219224078 I
 L := [0.4396034712, 1.545359967, -0.8928416048, -3.504544275]
     Lines := [-3.239253622 + 0.5 x, 0.7602488505 + 0.5 x, 8.699886068 + 0.5 x, 2.397905593 + 0.5 x]

          

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;