x0 = 0;
y0 = 0;
f = 0.245;
w = 0.6;
b = 0.15;
ndrive  = 200;
tstart = 0;
tend = ndrive*2*Pi/w ;

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

Plot[x[t] /. sol,  {t,  tend - 101,  tend - 1}];

tempx = x[tend] /. sol;
x0 = tempx[[1]]
tempy = y[tend] /. sol;
y0 = tempy[[1]]
tstart = tend;
w = 0.7;
tend = tend + ndrive*2*Pi/w;

ParametricPlot[{x[t], y[t]} /. sol, {t, tend - ndrive*2*Pi/w - 201, 
          tend - ndrive*2*Pi/w - 1}];

phase  =  
        0;

    px = Table[0, {i, 100, ndrive}];

    Do[{px[[i]] = 
            x[((i + 100) +   phase )*T ] /. sol[[1]]}, {i, 
          ndrive - 101}];

px;

py = Table[0, {i, 100, ndrive}];
Do[{py[[i]] = 
            y[((i + 100) +   phase )*T ] /. sol[[1]]}, {i, 
          ndrive - 101}];

py;

ListPlot[Transpose[{px, py}], PlotRange -> All, 
        PlotStyle -> PointSize[ .02]];