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

This is obviously a cooked up example.
The minimum is 0, but it is attained not only for x=y=z=1,  but also for  any x=1, y=a, z=a  (a>=1).