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 

@Markiyan Hirnyk 

A Google search gives

http://community.wolfram.com/groups/-/m/t/122095

@logiaco 

I would suggest to solve by hand this concrete example and show us how the result should be returned by Maple.
Note that you want only nontrivial solutions, so you may consider Determinant(A)=0.