(*Harmonic oscillator, exact solution of equations of motion*)
sol[ϕ0_,p0_,t_]={ϕ[t],p[t]}/.DSolve[{ϕ'[t]==p[t],p'[t]==-ϕ[t],ϕ[0]==ϕ0,p[0]==p0},{ϕ[t],p[t]},t][[1]]

(*Pendulum, numerical solution*)sol[ϕ0_,p0_]:=NDSolve[{ϕ'[t]==p[t],p'[t]==-Sin[ϕ[t]],ϕ[0]==ϕ0,p[0]==p0},{ϕ,p},{t,-10π,10π}][[1]]

(*Initial conditions oscillator: [1] points [2] region *)
points[1]={sol[0,0],sol[0,0.1],sol[0,0.2],sol[0,0.3],sol[0,0.4]};points[2]={sol[-.1,.4],sol[-.1,.5],sol[-.1,.6],sol[.0,.6],sol[.1,.6],sol[.1,.5],sol[.1,.4],sol[0,.4],sol[-.1,.4]};

(*plot range in units of π *) a=1/4;

pointplot[t_,set_]:=ListPlot[{Evaluate[{Mod[ϕ[t],2π,-π],p[t]}/.points[set]]},PlotMarkers->{Automatic,Large},AspectRatio->1,PlotRange->a{{-π,π},{-π,π}},AxesLabel->{"ϕ","p"},Joined->True];
Animate[pointplot[t,1],{t,0,2π},AnimationRunning->False,AppearanceElements->"ResetButton"]
(*change the argument <set> to switch between the different initial conditions*)

(*plot for longer times up till 8π to see the anharmonicity appear for the outer points*)

separatrix=Plot[{2Cos[t/2],-2Cos[t/2]},{t,-3.2,3.2},PlotRange->{-π,π},PlotStyle->{{Red,Thick},{Red,Thick}}]

piplot[s_]:=ParametricPlot[{{ s1 π,t π},{ s1 π,t π}}/.s1->s,{t,-1,1},PlotStyle->{Thick,Dashed,Orange}];(*Show[{piplot[1],piplot[-1]},PlotRange->{-π,π}]*)

stream=StreamPlot[{p,-Sin[ϕ]},{ϕ,-π,π},{p,-π,π},Axes->True, StreamScale->Medium,StreamStyle->Brown];

(*more initial conditions around the separatrix: [3] points [4] region *)
points[3]={sol[0,1.9],sol[0,1.95],sol[0,2.0],sol[0,2.05]};
points[4]={sol[-.1,1.9],sol[-.1,2.0],sol[-.1,2.1],sol[.0,2.1],sol[.1,2.1],sol[.1,2.0],sol[.1,1.9],sol[0,1.9],sol[-.1,1.9]};

(*change value for a to 1 to see a larger portion of phase space *) 
a=1;

Animate[Show[{stream,piplot[1],piplot[-1],separatrix,pointplot[tmax,3]},PlotRange->(a+0.1){{-π,π},{-π,π}}],{tmax,0,2π},AnimationRunning->False,AppearanceElements->"ResetButton"]