"Simulate N0 particles in a box"
SetDirectory[NotebookDirectory[]];
L=1.0;(*Box length*)
NP=10;(*number of particle*)
vref=0.1;(* reference speed *)

Ctimes=50;(*roughly the number of times the particle hitting the walls throughtout the simulation times*)timescale=2L/vref;
Deltat=0.01*timescale;
itlast=Ctimes*timescale/Deltat;
(*tlast=Deltat*itlast;*)


coordinates[1]={0,0};
coordinates[2]={0,L};
coordinates[3]={L,L};
coordinates[4]={L,0};
coordinates[5]={0,0};

"First, generate the box"
glbox=Graphics[Line[{Table[coordinates[n],{n,1,5}]}],AspectRatio->1]


(*Generate random initial conditions*)
Do[
v0[p]=vref;
v0x[p]=v0[p]*(2*Random[]-1);
v0y[p]=v0[p]*(2*Random[]-1);
Deltax[p]=Random[];
Deltay[p]=Random[];
{x[0,p] ,y[0,p]}={Deltax[p]*L,Deltay[p]*L};
,{p,1,NP}
];(*end do *)


lp=Graphics[
{PointSize[0.02],Point[Table[{x[0,p] ,y[0,p]},{p,1,NP}]]}
];

"Display initial configuration"
Show[{lp,glbox}]

For[np=1,np≤NP,np++,
xp[np]=x[0,np];
yp[np]=y[0,np]
];(*end for np*)

For[it=1,it≤itlast,it++,
For[np=1,np≤NP,np++,
xn[np]=xp[np]+v0x[np]*Deltat;(*the next x-position, xn, is calculated based on the previous postions xp*)
yn[np]=yp[np]+v0y[np]*Deltat;(*the next y-position, yn, is calculated based on the previous postions yp*)
x[it,np]=xn[np];
y[it,np]=yn[np];(*coordinates at time it are rocorded *)

(* choose boundary condition. Choose only one *)

(** Fixed boundary conditions **)
If[x[it,np]<0,x[it,np]=-x[it,np];v0x[np]=-v0x[np]];(*end if *)
If[y[it,np]<0,y[it,np]=-y[it,np];v0y[np]=-v0y[np]];(*end if *)
If[x[it,np]>L,x[it,np]=2L-x[it,np];v0x[np]=-v0x[np]];(*end if *)
If[y[it,np]>L,y[it,np]=2L-y[it,np];v0y[np]=-v0y[np]];(*end if *)
(** end of Fixed boundary conditions **)
(*
(** Periodic boundary conditions **)
If[x[it,np]<0,x[it,np]=L+x[it,np]];(*end if *)
If[y[it,np]<0,y[it,np]=L+y[it,np]];(*end if *)
If[x[it,np]>L,x[it,np]=x[it,np]-L];(*end if *)
If[y[it,np]>L,y[it,np]=y[it,np]-L];(*end if *)
(** end Periodic boundary conditions **)
*)

(* keep the record of xp, yp so that it can be used in the begining of the loop *)
xp[np]=x[it,np];yp[np]=y[it,np]
(* end of keep the record of xp, yp *)
](*end for np *)

];(*End  for it*)

"Output the simulationby use Graphics[Points] to visualise the simulation, which is a better option for such kind of simulation"

txy[it_]:=Table[{x[it,np],y[it,np]},{np,1,NP}];

GP[it_]:=Graphics[{PointSize[0.025],Point[txy[it]]},AspectRatio->1];

Manipulate[Show[GP[it],glbox],{it,0,itlast,1}]
(*Export["Nparticleinabox.flv",manip];*)