Motto: "In mathematics you don't understand things. You just get used to them" (John von Neumann).

A mathematician has an "abstract mind" so that in order to "see" R^13, he considers R^n and then takes n=13.
(This seems to be a joke, but basically it is not.)

Coming back more seriously to your example, observe first that in an abstract mathematical context such as Measure Theory,
Maple cannot help much. The visualization of a countable union of nested intervals should have been met in a Calculus course.

Anyway, in an abstract sigma-algebra, a visualization is out of the question; all we can do is to try special cases e.g. taking the space as finite, R, or R^2.

For the proof of the theorem 2.1 it is essentially enough to use the equivalence: 
sigma(A) = sigma(B)  iff  A subset sigma(B) and B subset sigma(A).

Visualizations are very seldom possible. Take for instance the proofs of the theorems:
1. Any Borel set is Lebesgue measurable.
2. There are subsets in R which are not Borel
(based on cardinality or on a Vitali-type construction).

restart;
JointNum := 3;                   
assign( seq(evaln(alpha[i]), i = 1 .. JointNum) = (0, 0, 0) );

If alpha is free (not assigned) then evaln can be omitted.

You cannot compute the indefinite integral numerically. Use:

seq(int(fN(x, 1.2), x=0..X), X=0..1,0.1);
0., 0.03008394009, 0.05955125696, 0.08778683199, 0.1141785882,   0.1381190848, 0.1590071982, 0.1762499145, 0.1892642546,  0.1974793487, 0.2003386758

The proof of the theorem: In a quadrilateral circumscribed to an ellipse, the line passing through the middle of the diagonals passes through the centre of the ellipse.
It is enough to have a circle instead of an ellipse (just take a projection). So, the proof:

restart;
x*cos(alpha[i])+y*sin(alpha[i])-1, x*cos(alpha[j])+y*sin(alpha[j])-1:
T:=unapply( eval([x,y], solve({%},[x,y])[]), [i,j]): # intersection for two tangents;
Slope:= P -> P[2]/P[1]:                              # the slope of the line OP
simplify( Slope((T(1,2)+T(3,4))/2) - Slope((T(1,3)+T(2,4))/2) );

        0

q.e.d.

msolve(5*x=7, 16);

                     {x = 11}

Yes, `evalf/EllipticCPi`(nu,k) is buggy. E.g.

evalf(EllipticCPi(0.985, 0.4));                 # crash
plot(EllipticCPi(n, 0.91), n=0.97..0.99); # discontinuous

Obvious workaround: use  EllipticPi(nu, sqrt(-k^2 + 1))
You may redefine

`evalf/EllipticCPi` := (nu,k) -> `evalf/EllipticPi`(nu, sqrt(1-k^2));

f:=CurveFitting:-Spline(my_array, x, degree=1):
b:=0.7: a:=solve(f=b,x):
linevert:=plot([a,y,y=0.1..0.9]):
display(my_graph, line, linevert);

restart;
NumPartStrict:=proc(n::integer[8], m::integer[8])::integer[8];
  option autocompile;
  local k;
  if m=0 then return `if`(n=0,1,0) fi;
  if m>=n then return 1 + thisproc(n,n-1) fi;
  add(thisproc(n-m+k-1,m-k),k=1..m)
end:

NumPartStrict(15,8);
                               13
NumPartStrict(125,75);
                            3179161
NumPartStrict(150,50);
                            14682366
 

You want the complex integral, not the line integral.

restart;
f:=exp(I*z)/z; Z:=r*exp(I*t):
J:=Int(eval(f, z=Z) * diff(Z,t), t=0..Pi);
value(J) assuming r>0;
evalf(eval(%,r=1/2)) = evalf(eval(J,r=1/2)); #check

For a non-polynomial expression we must define a custom order and use an inert version. For this example:

f:=x^2+x^y+1-z:
`%+`(sort([op(f)], (u,v) -> y in indets(u))[]);
#        x^y - z + 1 + x^2

restart;
ex:=int(diff(u(x,y,t),x)*v(x,y,t)+diff(v(x,y,t),x)*u(x,y,t),x): 
simplify(eval(ex, v=F/u)): eval(%, F=u*v);
#                     u(x, y, t) v(x, y, t)

It is always better to use explicit (or parametric) plots when possible.

plots:-animate(plot, [[sqrt(x^3/(2*a - x)), -sqrt(x^3/(2*a - x))], x = -5 .. 5, view=-5..5, color=red], a=-5..5);

If you really want implicitplot, just add the option signchange = false.

You could use inert Sum instead of sum and then, when needed do e.g.:
value(eval(a, N=4));