The procedure  InSphere  for any 4 points returns the center of the inscribed sphere and its radius. If the vertices of the tetrahedron are symbolic, then the result is also symbolic. Unfortunately, it takes quite a lot of time to do this calculation (see examples below). If the initial data (procedure arguments) are specified in float format, then the time is significantly reduced. You can play with this procedure for your own.

Download InSphere.mw

restart;
Eq1:=-T*sin(theta(t)) = m*(diff(X(t), t, t)+L*(diff(theta(t), t, t))*cos(theta(t))-L*(diff(theta(t), t))^2*sin(theta(t))):

Eq2:= T*cos(theta(t))-m*g = m*(diff(Y(t), t, t)+L*(diff(theta(t), t, t))*sin(theta(t))+L*(diff(theta(t), t))^2*cos(theta(t))):

numer(solve(Eq1, m)-solve(Eq2, m)):
simplify(-%/T)=0;
                     

When plotting replace  Im  with  'Im' .

plotmp_new.mw

No. For your case you need  unapply  command. An example:

a^2+b*c-c;
f:=unapply(%, a,b,c);
f(1,2,3);

plots:-display(plots:-implicitplot(y^4=y^2-x^2, x = 0.1 .. 0.3, y = 0.1 .. 0.3, color=blue), plot([[0.2, 0.2043096437]], style=point, symbol=solidcircle, color=red), view=[-0.5..0.5, -0.5..0.5]);
                 

Full view of the curve:

plots:-display(plots:-implicitplot(y^4=y^2-x^2, x = -1..1, y = -1..1, color=blue, gridrefine=3), plot([[0.2, 0.2043096437]], style=point, symbol=solidcircle, color=red), view=[-1..1, -1..1]);

Edit.
 

Since the probability density function is given explicitly, then all calculations are easy to make, simply by calculating the corresponding integrals without calling  Statistics  package. The code for the plottings can also be slightly simplified. 

Download NormaleVisualizationCourbe_new.mw

There are so many errors in the point (5) of your code that it is much easier to completely rewrite it. Of course the proposed version is particular. For example, it may happen that  M2[j, j]=0  or the system may be inconsistent or the equations may be dependent. For the general case, the code will be significantly more complicated.

Download LGS_Gauss-Jordan_new.mw

We can use  sum  instead of  add , then no problem:

sum(sum(1, i=1..j), j=1..4);    
S:=%sum(%sum(1, i=1..j), j=1..4);
value(S); 

Also in the general case:

sum(sum(1, i=1..j), j=1..n);
simplify(%);

 

Note that Maple automatically sorts the elements of your sets in the first list in ascending order. ListTools:-MakeUnique function removes duplicate elements from a list, keeping the order of the elements. Here is Ronan's example:

L1 := [{3, 5}, {4, 5}, {4, 8, 9}, {-7, 2, 3}]:
L2 := ListTools:-MakeUnique(op~(L1));

For your example, use the same commands. 

You should use 2D math Maple Notation (not LaTeX notation) for this. See  Help  ?plot,typesetting

You forgot to put a multiplication sign between 2 parentheses. In addition, it is not clear where  c+cos(x)  should stand, in the numerator or denominator. Maple has the same priority for division and multiplication. Therefore, for example, A/B*C  is understood as  (A*C)/B , and not as  A/(B*C) . Here are the calculations for both cases:

Expr1:=sin(x)/(a*b+a^2*sin(x)^2-d^2*cos(x)^2)*(c+cos(x));
Expr2:=sin(x)/((a*b+a^2*sin(x)^2-d^2*cos(x)^2)*(c+cos(x)));
int(Expr1, x);
int(Expr2, x);
                

Here the indefinite integrals are calculated. To calculate definite integrals, certain restrictions on the parameters are necessary, since a singularity can occur, i.e. zero in the denominator for your range of the integration.

A direct calculation (simbolic and numeric):

restart;
A:=simplify(int(ln(1+x)^3/(x+I), x = 0 .. 1));
B:=simplify(int(ln(1+x)^3/(x-I), x = 0 .. 1));
(A+B)/2;
simplify(evalf(%), zero);
 

The imagenary unit should be coded in Maple as  I  (not  i ).
Here is an example of solving the equation  f(w,B,V,k)=0  for  w  for specified values  B, V, k :

gm:=V->1/sqrt(1-V^2): T:=w-k*V: S:=w*V-k:
f:=unapply(I*B*gm(V)^3*T^3+gm(V)^2*T^2-I*gm(V)^3*T*S^2-1/3*I*B*gm(V)^3*S^3*T-1/3*gm(V)^2*S^2, w,B,V,k);
F:=(B,V,k)->fsolve(f(w,B,V,k), w, complex);
F(1, 0.5, 2);

If I understand correctly, you have an equation with several parameters and you want to plot a graph that shows how the roots change depending on the change of parameters. To do this, you need to solve the equation numerically using  fsolve  command. Since the roots are not one, but several, then the number of plots should be the same as the number of roots. In the example below, we explore how the imaginary part of the roots of a quadratic equation changes depending on the change in parameters  p  and  q :

f:=(x,p,q)->x^2+p*x+q:
F:=(p,q)->fsolve(f(x,p,q), complex);
P1:=plot3d(Im('F'(p,q)[1]), p=-10..10, q=-10..10, color=grey, numpoints=10000, axes=normal):
P2:=plot3d(Im('F'(p,q)[2]), p=-10..10, q=-10..10, color=grey, numpoints=10000, axes=normal):
plots:-display(<P1 | P2>);

The graphs clearly show that the region where the imaginary part is 0, i.e. the roots are real, separated by the parabola  p^2-4*q=0  (i.e. the discriminant is 0), and the values of the complex roots are conjugate.

Edit.

Download maple17_new.mw