@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

@acer 

The plot component has a delay option (default 100 ms) which seems to be ignored. Do you know why?

Does not help. I wonder if it's only my configuration.

Edit. You have posted probably an older version where you forgot to change
cstep to step.
After the correction it works.

@Carl Love 

@Carl Love 

The code does not work in Maple (2016, 64 bit, Win7).

@Carl Love 

Carl, on my computer the animation does not work, all the frames are identical.

PS. How have you included the gif file? In my case a single frame was copied this time.

@Carl Love 

However I think that such an animation could be possible, updating continuously a plot component (the compiled step seems to be fast enough for this). I let somone else to try.

Edit. acer already found a simpler solution using Explore.

BTW, the Explore command is now very powerful and probably should be used more often.

@Carl Love 

You are probably right here, but in our situation the order of the enties is not important provided that it is coherent (and it should be).

@Carl Love 

Iterator uses a one-dimensional Array, for which entries is documented. But it works indeed for Vectors too.

@taro 

Just replace beta[1,2]  with beta__12 or other symbol.

Smith form can be used to find ker(A) in all cases.
An example:

A := Matrix([[1,2*p,2*p^2+2*p],[1,6*p,6*p^2+6*p],[1,2,p]]):
S := SmithForm(A); 

(U,V):=SmithForm(A, p, output=['U', 'V']):   #U.A.V=S
# if p^2+2*p<>0 then X=0
# if p^2+2*p=0, ker(A.V) is generated by
K:=<0,0,1>;
# so, ker(A) is generated by
X0:=V.K;


# Check:
A.X0;  #  = 0  (since p^2+2*p=0)

@logiaco