(* Rownanie wahadla tlumionego z wymuszeniem periodycznym *)
(* alpha - kat, a - wspolczynnik tlumienia, w0 - czestosc wlasna, A - amplituda sily wymuszajacej, w1 - czestosc sily wymuszajacej *)
eq = alpha''[t] + a alpha'[t]+w0^2 Sin[alpha[t]] == A Cos[w1 t];


(* Parametry *)
a = 1/10; w0 = 1; w1 = 1; A = 1;

(* Czas do ktorego badany jest ruch wahadla *)
tmax=200;

(* Liczba zwiazana z liczba nawiniec (0 jesli nie mozna rozpoznac liczby nawiniec) *)
atraktor[al_,alp_]:=Module[{sol=NDSolve[{eq,alpha[0]==al, alpha'[0]==alp},alpha,{t,0,tmax}],tbl},
tbl=Quotient[Table[alpha[100+2 Pi i]/.sol,{i,15}],Pi];
If[Take[tbl,-5]==Table[Last[tbl],{5}],Last[tbl][[1]],0]
];

(* Badanie fragmentu plaszczyzny danych poczatkowych z wykorzystaniem wszystkich dostepnych rdzeni *)
DistributeDefinitions[eq,a,w0,w1,A,tmax,atraktor];
SetSharedVariable[y];
For[y=0,y<=2,y=y+0.0005,row=ParallelTable[{x,y,atraktor[x, y]},{x,0,2,0.0005}];Do[Print[row[[i,1]],"\t",row[[i,2]],"\t",row[[i,3]]],{i,1,Length@row}]]
Exit[]