x0 = 0;
y0 = 0;
f = 0.26;
T = 2*Pi;
b = 0.15;
ndrive  = 1000;
tstart = 0;
tend = ndrive*T ;

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

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

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

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]];