Clear["Global`*"]

n = 64; (*number of non-endpoint oscillators*)
α = 1; (*nonlinearity parameter*)
tMax = 2000; (*maximum time*)
f0 = 3;(*all energy goes into the "f0-th" mode*)

Clear[f]
f[y_,α_] = y + α y^2; (*f is the derivative of the potential function*)
Clear[x]
s = Join[{0},Table[x[i,t],{i,1,n-1}],{0}];
st = Table[x[i,t],{i,1,n-1}];
Clear[sys]
sys = Table[D[s[[i+1]],{t,2}]==f[s[[i+2]]-s[[i+1]],α]-f[s[[i+1]]-s[[i]],α],{i,1,n-1}];
sysInit = Table[Evaluate[(s[[i+1]])/.{t->0}] ==Sin[((f0) π)(i)/n],{i,1,n-1}];
sysInitV=Table[Evaluate[(D[s[[i+1]],t])/.{t->0}]==0,{i,1,n-1}];
Clear[sol]
sol = NDSolve[Join[sys,sysInit,sysInitV],st,{t,0,tMax},MaxSteps->100000];

movieFPU=Animate[ListPlot[Join[{0},Flatten[Evaluate[Table[x[i,t],{i,1,n-1}]/.sol]/.t->r],{0}],PlotRange->{{1,n+1},{-1,1}},Joined->False,PlotStyle->{PointSize[.011],RGBColor[.5,.8,1]},Background->Black,ImageSize->700],{r,0,tMax,1},AnimationRate->17]