How to plot poincare map as shown in below image for Duffing equation?

Question

I am trying to plot phase diagram and poincare map. But I cannot get the poincare map as shown in the image below

Phase Space

sol = NDSolve[{v'[t] == 
     0.320 x[t] - 1.65 x[t]^3 - 0.005*v[t] + 0.855 Cos[1.2*t], 
    x'[t] == v[t], x[0] == 0, v[0] == 0}, {x, v}, {t, 0, 1500}];
ParametricPlot[{x[t], v[t]} /. sol, {t, 200, 1000}, 
 AxesLabel -> {"x", "v"}, PlotRange -> Full, PlotStyle -> LightGray, 
 Axes -> False, Frame -> True, 
 FrameTicksStyle -> Directive[Black, 20], ImageSize -> {700, 350}, 
 AspectRatio -> Full]

Poincare Map

poincare[A_, gamma_, omega_, ndrop_, nplot_, 
  psize_] := (T = 2*Pi/omega;
  g[{xold_, vold_}] := {x[T], v[T]} /. 
    NDSolve[{v'[t] == 
        0.320 x[t] - 1.65 x[t]^3 - gamma*v[t] + A*Cos[omega*t], 
       x'[t] == v[t], x[0] == xold, v[0] == vold}, {x, v}, {t, 0, 
       T}][[1]];
  lp = ListPlot[Drop[NestList[g, {0, 0}, nplot + ndrop], ndrop], 
    PlotStyle -> {PointSize[psize], Black}, Axes -> False, 
    Frame -> True, FrameTicksStyle -> Directive[Black, 20], 
    PlotRange -> All, AxesLabel -> {"x", "v"}, 
    ImageSize -> {700, 350}, AspectRatio -> Full])

poincare[0.855, 0.005, 1.2, 1000, 200, 0.01]

I want diagram as shown in below image(The phase diagram will be different for my code)

Answer 1

If you run your first code for the same amount of time as the Poincare section, you'll see that they do agree (and that they're not particularly interesting parameter values):

T = 2*Pi/1.2;

sol = NDSolve[{v'[t] == 0.320 x[t] - 1.65 x[t]^3 - 0.005*v[t] + 0.855 Cos[1.2*t], 
  x'[t] == v[t], x[0] == 0, v[0] == 0}, {x, v}, {t, 0, 1200 T}][[1]];

pp = ParametricPlot[{x[t], v[t]} /. sol, {t, 1000 T, 1200 T}, 
  AxesLabel -> {"x", "v"}, PlotRange -> Full, PlotStyle -> LightGray, 
  Axes -> False, Frame -> True, 
  FrameTicksStyle -> Directive[Black, 20], ImageSize -> {700, 350}, 
  AspectRatio -> Full, MaxRecursion -> 7];

poincare[0.855, 0.005, 1.2, 1000, 200, 0.005];
Show[pp, lp]

Answer 2

Not an answer, but an essential warning !

For t->1500 you have more than 300 oscillations of x[t]. Even at best Workingprecision numerical NDSolve will give total wrong result for higher t. See examples. A workaround could be solutions with DSolve (but my MMA version doesn't find any.)

sol1 = NDSolve[{v'[t] == 
 0.320 x[t] - 1.65 x[t]^3 - 0.005*v[t] + 0.855    Cos[1.2*t], 
x'[t] == v[t], x[0] == 0, v[0] == 0}, {x, v}, {t, 0, 1500}, 
MaxSteps -> 10^5];

sol2 = NDSolve[{v'[t] == 
  0.320 x[t] - 1.65 x[t]^3 - 0.005*v[t] + 0.855 Cos[1.2*t], 
 x'[t] == v[t], x[0] == 0, v[0] == 0} // Rationalize[#, 0] &, {x, 
v}, {t, 0, 1500}, MaxSteps -> 10^5, WorkingPrecision -> 25];

sol3 = NDSolve[{v'[t] == 
 0.320 x[t] - 1.65 x[t]^3 - 0.005*v[t] + 0.855 Cos[1.2*t], 
x'[t] == v[t], x[0] == 0, v[0] == 0} // Rationalize[#, 0] &, {x, 
v}, {t, 0, 1500}, MaxSteps -> 10^6, WorkingPrecision -> 50]

{xsol1, vsol1} = {x, v} /. First@sol1;
{xsol2, vsol2} = {x, v} /. First@sol2;
{xsol3, vsol3} = {x, v} /. First@sol3;

Plot[{xsol1[t], xsol2[t], xsol3[t]}, {t, 300, 400}, 
PlotStyle -> {Blue, Red, Green}]