(*The Duffing Equation:example of hysteresis 
(damped,driven oscillator with force f[x]=-k x+eps x^3)*)

(*enter initial conditions *)

x0=0;(*initial x-value *)
y0=0;(*initial y-value *)
f=0.01;(*magnitude of driving force *)
w=1;(*driving frequency *)
gamm=0.15;(*damping constant *)
ndrive=1000;(*number of driving periods *)
tstart=0;(*initial time *)
tend=ndrive*2*Pi/w;   (*final time *)

(* Set up to solve the Duffing equation
rewrite second order d.e.as two first order d.e.s  *)

(*The Duffing equation:
x''+w0^2 x-w0^2 x^3/lambda^2+b/m x'=F0/m*cos[w t] 
(after setting m=1,lambda=1,w0^2=1) reduces to:
x'=y and 
y'=-gamm y+x-x^3+f cos[w t] where gamm=b/m   *)

sol=NDSolve[{x'[t]==y[t],
           y'[t]==-gamm*y[t]+x[t]-(x[t])^3+f*Cos[2*w*t],
           y[0]==y0,x[0]==x0},
             {y,x},{t,tstart,tend},MaxSteps->200000]

(*Plot x[t] well after any transient motion, i.e. after 1100 driving periods *)

Plot[x[t]/.sol,{t,1100,1150}]

(*  Create the phase space plot of {x,x'(=y)}    *)

ParametricPlot[{x[t],y[t]}/.sol,{t,1100,1600}]

(* Poincare Section *)

phase = 0;  (*phase shift *)
T = 2 Pi/w;  (*period of driving force *)
px=Table[0,{i,100,ndrive}]; (*initialise px with zero entries *)
(*extract p[x] fpr same place in every cycle of the driving period *)
Do[{px[[i]]=x[((i+100)+  phase )*T ]/.sol[[1]]},{i,ndrive-101}];

py=Table[0,{i,100,ndrive}];(*initialise py with zero entries *)
(*extract p[x] fpr same place in every cycle of the driving period *)Do[{py[[i]]=y[((i+100)+  phase )*T ]/.sol[[1]]},{i,ndrive-101}];

(* Plot Poincare section, i.e. phase plot at same point in time at every driving period *)

ListPlot[Transpose[{px,py}],PlotRange->All,PlotStyle->PointSize[0.02]]