@Carl Love 

Of course, thank you.

@Mariusz Iwaniuk 

Try

plots[odeplot](de, [x, diff(y(x),x,x)^2 - 1], x=0..1);

OK,  so, if rating(X)=x, rating(Y)=y  ==>

Prob(X_beats_Y) = 1 / (1 + 10^((y-x)/400))

@Christopher2222 

Do you know how the ELO ratings are computed?
If two teams X,Y have the ratings x,y what is the probability for X to win?

(I see in the worksheet that Gumbel distribution is used, but I wonder why it was chosen.)

@Lola 

Your function mi() is discontinuous in (tau,T) due to charfcn.This will creat problems for NLP. Also, charfcn is probably not the best choice for numerics.
I would recommend to investigate theoretically the problem.

For vizualization I suggest also e.g.:

Explore(plot(mi(10, 1, 5, 2, 2, 6, T, tau, 0.2, 0.2, 0.5),T=0..20,view=-1..1),tau=0.0..20.0,animate);

@Markiyan Hirnyk 

I don't know of course which summation method was used. But you can easily verify whether the result agrees with the one given by Euler's summation method.

@Zeineb 

That's simply because your new series is divergent.
Try e.g.

add(evalf(g(n)/n!),n=0..500);

@Markiyan Hirnyk 

What is aside? The Fibonacci sequence is a linear combination of geometric sequences, isn't it? Have you asked something else?

The geometric series Sum(z^n,n=0..infinity) converges to 1/(1-z) in Euler's sense,

see https://en.wikipedia.org/wiki/Euler_summation

@Markiyan Hirnyk 

@Markiyan Hirnyk 

But note that in this case it is actually used
Sum(x^n, n=0..infinity) = 1/(1-x);   # for |x|>1   !!

@Carl Love 

Here is an "inline" example as I understand it; I am not sure whether this is what acer suggested.

n:=10:
XY:=Matrix(n,2, (i,j)->`if`(j=1,i,i^2),datatype=float[8]):

updatexy:=proc(XY::Matrix(datatype=float[8]))
option autocompile;
local i::integer[4];
for i to n do XY[i,2]:=XY[i,2]*0.999 od;
end:

P:=proc(XY::Matrix(datatype=float[8]))
updatexy(XY);
PLOT(CURVES('XY', COLOUR(RGB, .47058824, 0., 0.54901961e-1, _ATTRIBUTE("source" = "mathdefault"))), STYLE(POINT), AXESLABELS("", ""), VIEW(1 .. n, 1 .. n^2)); 
 # se pare ca nu conteaza 'XY'
end:

P1:=proc(XY::Matrix(datatype=float[8]))
updatexy(XY);
plot(XY, style=point, view=[1..n,1..n^2]);
end:

t:=time[real]():
Explore(P(XY),parameters=[[k=1..10,shown=false,animate]],numframes=1000,autorun);
time[real]()-t; fps=1000/%;

I am not sure whether the XY matrix is passed to the plot engine by reference or it is "expanded".
Anyway, for n=10000 the animation is considerably slower even if the computation time for  the compiled procedure updatexy is rather small (if updatexy is set to do nothing, the animation speed is almost the same).
My guess is that plotting a large XY matrix is slow even if the plot structure is updated inline.
On my computer for n=10000 Explore can dispay approx 17 frames / second. For n=10 it is 62 fps.
For P1 and n=10000 ==> 8.6 fps (twice slower, because plot generates first a PLOT structure)

I suppose that there is the possibility that some dancers escape from the dance floor; maybe even get lost in the space :-). The code could help finding them!

@Carl Love 

(Carl's idea with acer's excellent method using Explore). Of course it could be enhanced with other params ...

n:=4000: rechoose:=10:
x,y,z:='LinearAlgebra:-RandomVector(n,generator=-1.0 ..  1.0,datatype=float[8])' $3:
FE:=proc() global F,E; #friend,enemy
  F,E:='LinearAlgebra:-RandomVector(n,generator=1..n,datatype=integer[4])' $2;
end:
FE():

step:=proc(x::Vector(datatype=float[8]),y::Vector(datatype=float[8]),z::Vector(datatype=float[8]),
      E::Vector(datatype=integer[4]),F::Vector(datatype=integer[4]),mu::integer[4])
option autocompile;
local i::integer[4],ex::float[8],ey::float[8],ez::float[8],ed::float[8],
                    fx::float[8],fy::float[8],fz::float[8],fd::float[8];
to mu do
for i to n do
  ex:=x[E[i]]-x[i];ey:=y[E[i]]-y[i]; ez:=z[E[i]]-z[i]; ed:=sqrt(ex^2+ey^2+ez^2)+0.01;
  fx:=x[F[i]]-x[i];fy:=y[F[i]]-y[i]; fz:=z[F[i]]-z[i]; fd:=sqrt(fx^2+fy^2+fz^2)+0.01;
  x[i] := 0.995*x[i] - 0.01*ex/ed + 0.02*fx/fd;
  y[i] := 0.995*y[i] - 0.01*ey/ed + 0.02*fy/fd;
  z[i] := 0.995*z[i] - 0.01*ez/ed + 0.02*fz/fd;
od; od;
end:

P:=proc(mu)
  step(x,y,z,E,F,mu);
  if rand(1..100)()<=rechoose then FE() fi;
  plots:-pointplot3d(x,y,z,axes=boxed,view=[-2..2,-2..2,-2..2],
       color="blue",symbolsize=3,tickmarks=[0,0,0]); #symbol=circle
end proc:

Explore(P(mu),parameters=[[xx=0..1,shown=false,animate], [mu=5..20,animate=false]],
        initialvalues=[mu=10],autorun,loop,numframes=100);

Probably "true or" just remained there after some test and should be removed.

xx is a dummy parameter; Explore needs a at least 1 animated parameter in order to produce animation. But xx could be removed in the Starlings(...) call.

@Carl Love 

Thank you for the detailed and competent answer.
(I have already tried Sleep and it works as expected.)
@acer