Questa mappa genera le soluzioni

solT[{x0_,v0_}]:=Module[{},
{x[t],v[t]}/.NDSolve[{x'[t]==v[t],v'[t]==Cos[t]-0.1 v[t]-Sin[x[t]],x[0]==x0,v[0]==v0},{x[t],v[t]},{t,0,60π}][[1]]]/.t->60π

Questa è la tabella dei dati iniziali

array=Table[{i,j},{j,-2.,2.,0.0052083333333334},{i,-N[3π],N[2π],0.0153398078788565}];

Questa è la tabella delle evolute dopo 60 π

{tempo,arrayT}=Map[solT,array,{2}]//Timing;

tempo

Definiamo una funzione che associa i valori di colore rosso, verde e blu ad i punti che tendono verso le "prime" 3 orbite periodiche. Bianco a quelli che non si riesce a decidere

col[{x_,v_}]:=Which[
(((x+2.046291506120066-2π))^2)+(((v-0.39269675159911344))^2)<0.1,0,
(((x+2.046291506120066))^2)+(((v-0.39269675159911344))^2)<0.1,1,
(((x+2.046291506120066+2π))^2)+(((v-0.39269675159911344))^2)<0.1,2,
1==1,3]

Facciamo il disegno

ArrayPlot[Map[col,arrayT,{2}],
ColorRules->{0->Red,1->Green,2->Blue,3->White}]