Use the Student:-VectorCalculus package and the LineInt command for this:

restart;
with(Student:-VectorCalculus):
LineInt(VectorField(<x^2-y^2+x, -2*x*y-y>), Path(<t^2, t>, t = 0 .. 1));
LineInt(VectorField(<x^2-y^2+x, -2*x*y-y>), Path(<t^2, t>, t = 0 .. 1),output=plot,fieldoptions=[fieldstrength = fixed]);

The expression for your function can be greatly simplified to the fraction  (x-4)/x  by reducing the numerator and denominator by  x^2-9 . But at the same time you lose information that the original expression in points  x=-3  and  x=3  will be  0/0  (these will be the points of removable discontinuity) .
For better plotting, it is useful to increase the range on the x-axis, it is reasonable to limit the range on the y-axis. Additionally, the horizontal asymptote, vertical asymptote and points of removable discontinuity are shown in blue on the plot below:

restart;
f:=x->(x^3 - 4*x^2 - 9*x + 36)/(x^3 - 9*x);
g:=unapply(simplify(f(x)),x);
Lines:=plot([f(x), [0,t, t=-10..10],1], x=-7..7, y=-6..7, linestyle=[1,3,3], thickness=[2,0,0], color = [red,blue$2], discont ):
Points:=plot([[-3,g(-3)],[3,g(3)]], style=point, symbol=circle, color=blue, symbolsize=14):
plots:-display(Lines,Points);

Addition. If we write  (x - 4)/x = 1 - 4/x , then we see that the graph of your function will be symmetrical about the point  (0, 1)  (this graph is called a hyperbola).

a. To plot a semicircle graph, it is better to use parametric equations, which make it easy to plot any part of the circle.
b. Obviously, the rotation results in a ball of radius 2 centered at the origin.
c. To calculate the volume of this ball, you can use the explicit equation of the given semicircle and the penultimate formula from your list.
 

restart;
plot([2*cos(t),2*sin(t),t=0..Pi], scaling=constrained);
f(x):=sqrt(4-x^2):
Pi*int(f(x)^2, x=-2..2);

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]);