SetPendParam[driveforce_, drivefreq_, dampen_, phase_] := 
        Module[{q = driveforce, p = phase, b = dampen, 
            w = drivefreq}, 
f1[t_, x_, y_, z_] = 
            y; 
f2[t_, x_, y_, z_] = 
            q * Cos[z + p] - Sin[x] - b * y; 
f3[t_, 
              x_, y_, z_] = 
            w; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], f3[t, x, y, z]};];

Runge[initt_, finalt_, initx_, inity_, 
          steps_] := 
Module[{it = initt, ft = finalt, j, 
            m = steps}, 
h = (ft - it)/
              m; 
T = 
            Table[a + (j - 1) h, {j, 1, 
                m + 1}]; 
R = {initx, inity, 
              initt}; 
P = 
            Table[{}, {m}]; 
PP = 
            Table[{}, {m}]; 
For[j = 1, 
            j < m + 1, j++, 
k1 = 
              h F[T[[j]], R]; 
k2 = 
              h F[T[[j]] + h/2, R + k1/2]; 
k3 = 
              h F[T[[j]] + h/2, R + k2/2]; 
k4 = 
              h F[T[[j]] + h, R + k3]; 
Rnew = 
              R + 1/6 (k1 + 2 k2 + 2 k3 + k4); 
R = 
              Rnew;
]; 
Return[
            R];
];

other  =  Table[{0, 0}, {1}];
Do[
{α = 
            Random[Real, {0.8, 
                1.8}]; 
SetPendParam[α, 2/3, 0.5, 
            0]; 
a = 
            Random[Real, {0, 2}]; 
aa = 
            Random[Real, {0, 2}]; 
Runge[0, 125, a, aa, 
            2500]; 
β = 
            Mod[R[[1]], 
                2  π] - π; 
PPP  = {{α, β}}; 
other  =  
            Join[PPP, 
              other];
}
, {30000}];
ListPlot[other, PlotStyle -> PointSize[0.001], 
        PlotRange -> {{0.8, 
              1.8}, {-π, π}}];
other  =  Table[{0, 0}, {1}];
Do[
{α = 
            Random[Real, {0.8, 
                1.8}]; 
SetPendParam[α, 2/3, 0.5, 
            0]; 
a = 
            Random[Real, {0, 2}]; 
aa = 
            Random[Real, {0, 2}]; 
Runge[0, 125, a, aa, 
            1250]; 
β = 
            R[[2]]; 
PPP  = {{α, β}}; 
other  =  
            Join[PPP, 
              other];
}
, {30000}];