Clear[wd,g,q,w0,th0,ndrive,tlim,sol,th,w,t]

Driven Damped Pendulum

wd = 0.8;                               (* Drive freqency          *)
f = 0.3;                                  (* Drive force amplitude   *)
c= 0.15;                                 (* Damping factor          *)
w0 =0;                                      (* Initial velocity        *)
th0 =0;                                  (* Initial angle           *)
ndrive = 400;                         (* Number of drive cycles  *)
tlim = ndrive*2 Pi/wd      (* Total time              *)

(* Set up to solve the driven, damped pendulum, M&T eqn 4.43 *)

sol = NDSolve[{th'[t] == w[t],
		                              w'[t]  == -c * w[t] -Sin[th[t]] + f Cos[wd t],
		                             w[0]  == w0,
		                             th[0] == th0},
	                               {w,th},{t,tlim}, MaxSteps ->610000];

(* Plot the pendulum angle as a function of time *)

Plot[th[t]/.sol,{t,1100,1200}]

(* Plot the pendulum's angular speed a function of time *)

Plot[w[t]/.sol,{t,1100,1200}]

(* Plot the phase diagram for the pendulum *)

ParametricPlot[{th[t],w[t]}/.sol,{t,800,1000}]

(* rescale the phase diagram so that the angle lies between -pi and +pi. *)

ParametricPlot[{Mod[th[t]+π,2 π]-π,w[t]}/.sol,{t,800,1200}]

(* Poincare Section *)
(* Plotting the Poincare section by looking at the phase diagram at the same part 
of the driving period each cycle.  *)
phase = 0;
pw = Table[0,{i,20,ndrive}]; (*initialise the pw table *)
Do[{pw[[i]] = w[((i+20)+  phase )*2 Pi/wd ]/.sol[[1]]},{i,ndrive-21}];

pth = Table[0,{i,20,ndrive}]; (*initialise the pth table *)
Do[{pth[[i]] = Mod[th[((i+20) + phase )*2 Pi/wd ]+ Pi,2 Pi] - Pi/.sol[[1]]},{i,1,ndrive-21}]

ListPlot[Transpose[{pth,pw}],PlotRange->All,PlotStyle->PointSize[0.02]]