(*Simulation Parameters*)
NumParticles=500;
Lmin=-1;
Lmax=1;
Pmin=-.01;
Pmax=.01;
WallXmin=Lmin;
WallXmax=0;
WallPmin=Pmin;
WallPmax=Pmax;
m=1;
dt=.1;
tmax=2000;

(*Simulation Code - Don't Change*)
t=Range[0,tmax,dt];
length=Length[t];

Do[
XPPoints[i]=Module[{t},Table[If[t<1,
{X[0]=RandomReal[{WallXmin,WallXmax}],P[0]=RandomReal[{WallPmin,WallPmax}]},
If[Abs[X[t-1]]+P[t-1]/m*dt≥Lmax,{X[t]=2*Lmax*P[t-1]/Abs[P[t-1]]-X[t-1]-P[t-1]/m*dt,P[t]=-P[t-1]},{X[t]=X[t-1]+P[t-1]/m*dt,P[t]=P[t-1]}]],
{t,0,length}]],
{i,1,NumParticles}
]

data=Table[XPPoints[n],{n,NumParticles}]  //Transpose;

Animate[ListPlot[List/@data[[t]],Axes->True,AxesLabel->{X,P},PlotLabel->"Free Particle Phase Space", PlotRange->{{Lmin,Lmax},{Pmin,Pmax}}],{t,1,Length[data],1},AnimationRate->500,AnimationRunning->True]