SetPendParam[driveforce_, drivefreq_, dampen_, phase_] := 
        Module[{q = driveforce, b = dampen, p = phase, 
            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_, initx_, inity_, stepsper_, 
          numb_] := 
Module[{it = initt, j, 
            stepper = stepsper, num = numb}, 
period = 
            2.0*N[π, 20]*(3/2); 
h = 
            (period/stepper); 
third  =  
            Table[{}, {stepper}]; 
thirdb  =  
            Table[{}, {num}]; 
T = 
            Table[(it + (j - 1)* h), {j, 1, 
                stepper}]; 
R = {initx, inity, 
              it}; 
Do[{For[j = 1, 
                  j < stepper + 1, j++, 
k1 = 
                    h*F[j*T[[j]], R]; 
k2 = 
                    h*F[j*(T[[j]] + h/2), 
                        R + k1/2]; 
k3 = 
                    h*F[j*(T[[j]] + h/2), 
                        R + k2/2]; 
k4 = 
                    h*F[j*(T[[j]] + h), 
                        R + k3]; 
Rnew = 
                    R + 1/6 (k1 + 2 k2 + 2 k3 + 
                            k4); 
R = 
                    Rnew; 
second = 
                    R; 
third[[j]] = 
                    second
]; 
thirdb[[k]] = third
}, {k, 1, 
                num}]
Do[
              thirdb[[g + 1]] = 
                Join[thirdb[[g]], thirdb[[g + 1]]], {g, 1, 
                num - 1}]; 
thirdb  =  Drop[thirdb, 
                num - 1][[1]];
];

PlotTheta[fvalue_, bvalue_, steps_, 
        rate_] := 
Module[{fff  =  fvalue, 
          bbb  =  bvalue,  ssteps  =  steps,  
          rrate  =  rate}, 
SetPendParam[fff, 2/3, bbb, 
          0]; 
fourth  =  
          Table[{}, {steps * rate}]; 
Runge[0, 1, 1, 
          rate, steps]; 
Do[
          fourth[[k]]  =  {thirdb[[k]][[3]], thirdb[[k]][[1]]}, {k, 1, steps * rate}]; 
ListPlot[
          fourth];
]

CalculateTheta[fvalue_, bvalue_, steps_, rate_, initx_, 
        inity_] := 
Module[{fff  =  fvalue, 
          bbb  =  bvalue,  ssteps  =  steps,  rrate  =  rate, 
          initth  =  initx,  
          initw  =  inity}, 
SetPendParam[fff, 2/3, bbb, 
          0]; 
fourth  =  
          Table[{}, {steps * rate}]; 
Runge[0, initth, 
          initw, rate, steps]; 
Do[
          fourth[[k]]  =  thirdb[[k]][[1]], {k, 1, 
            steps * rate}]; 
Return[
          fourth];
]

CompareTheta[fvalue_, bvalue_, steps_, rate_, initx_, inity_, 
        delta_] := 
Module[{fff  =  fvalue, 
          bbb  =  bvalue,  ssteps  =  steps,  rrate  =  rate, 
          initth  =  initx,  
          initw  =  inity, δ = delta}, 
blah  =  
          CalculateTheta[fff, bbb, ssteps, rrate, initth, 
            initw]; 
blah2  =  
          CalculateTheta[fff, bbb, ssteps, rrate, initth + δ, 
            initw + δ]; 
blah3  =  
          Table[{0}, { .9*ssteps*rrate}]; 
Do[
          blah3[[i]]  =  
            Log[Abs[blah2[[i +  .1*ssteps*rrate]] - 
                  blah[[i +  .1*ssteps*rrate]]]], {i, 
            1,  .9*ssteps*rrate}]; 
gp  = 
          ListPlot[blah3, 
            DisplayFunction -> Identity];
]

CalculateOmega[fvalue_, bvalue_, steps_, rate_, initx_, 
        inity_] := 
Module[{fff  =  fvalue, 
          bbb  =  bvalue,  ssteps  =  steps,  rrate  =  rate, 
          initth  =  initx,  
          initw  =  inity}, 
SetPendParam[fff, 2/3, bbb, 
          0]; 
fourth  =  
          Table[{}, {steps * rate}]; 
Runge[0, initth, 
          initw, rate, steps]; 
Do[
          fourth[[k]]  =  thirdb[[k]][[2]], {k, 1, 
            steps * rate}]; 
Return[
          fourth];
]

CompareOmega[fvalue_, bvalue_, steps_, rate_, initx_, inity_, 
        delta_] := 
Module[{fff  =  fvalue, 
          bbb  =  bvalue,  ssteps  =  steps,  rrate  =  rate, 
          initth  =  initx,  
          initw  =  inity, δ = delta}, 
blah  =  
          CalculateOmega[fff, bbb, ssteps, rrate, initth, 
            initw]; 
blah2  =  
          CalculateOmega[fff, bbb, ssteps, rrate, initth + δ, 
            initw + δ]; 
blah3w  =  
          Table[{0}, { .9*ssteps*rrate}]; 
Do[
          blah3w[[i]]  =  
            Log[Abs[blah2[[i +  .1*ssteps*rrate]] - 
                  blah[[i +  .1*ssteps*rrate]]]], {i, 
            1,  .9*ssteps*rrate}]; 
gp  = 
          ListPlot[blah3w, 
            DisplayFunction -> Identity];
]

CompareTheta[1.5, 0.5, 10, 100, 1, 1, 
        0.01];
Show[gp, 
        DisplayFunction -> $DisplayFunction];
CompareOmega[1.5, 0.5, 10, 100, 1, 1, 
        0.01];
Show[gp, DisplayFunction -> $DisplayFunction];

Lyapunov[fvalue_,  bvalue_,  cycles_,  rate_,  delta_,  
        number_] := 
Module[{fff  =  fvalue,  
          bbb  =  bvalue,  ccycles  =  cycles,  
          rrate  =  rate,  δ =  delta,  
          numb  =  number}, 
lyapunov1  =  
          Table[0, {numb}]; 
lyapunov2  =  
          Table[0, {numb}]; 
phasediff1  =  
          Table[0, {rrate}]; 
phasediff2  =  
          Table[0, {rrate}]; 
Do[{a = 
              Random[Real, {0, 1}]; 
aa = 
              Random[Real, {0, 1}]; 
CompareTheta[fff, bbb, 
              ccycles, rrate, a, 
              aa, δ]; 
CompareOmega[fff, bbb, 
              ccycles, rrate, a, aa, δ]; 
Do[
              phasediff1[[j]]  =  blah3[[j]], {j, 1, 
                rrate}]; 
Do[
              phasediff2[[j]]  =  blah3w[[j]], {j, 1, 
                rrate}]; 
lyapunov1[[k]]  =  
              D[Fit[phasediff1, {1, x}, x], 
                x]; 
lyapunov2[[k]]  =  
              D[Fit[phasediff2, {1, x}, x], x];
}, {k, 1, 
            numb}]
]

Lyapunov[1.5, 0.5, 20, 50,  .1, 2]
lyapunov1
lyapunov2

PlotLyapunov[bvalue_, cycles_,  rate_,  delta_,  number_, 
        tot_] := 
Module[{bbb  =  bvalue,  
          ccycles  =  cycles,  rrate  =  rate,  δ =  delta,  
          numb  =  number, total = tot}, 
other1  =  
          Table[{0, 0}, {1}]; 
other2  =  
          Table[{0, 0}, {1}]; 
Do[{α = 
              Random[Real, {0.9, 1.1}]; 
Lyapunov[α, 
              bbb, ccycles, rrate, δ, 
              numb]; 
β1  =  
              Mean[lyapunov1]; 
β2  =  
              Mean[lyapunov2]; 
P  =  {{α, β1}}; 
PP  =  {{α, β2}}; 
other1  =  
              Join[P, other1]; 
other2  =  
              Join[PP, 
                other2];
}, {total}]; 
other1  =  
          Drop[other1,  -1]; 
other2  =  
          Drop[other2,  -1]; 
ListPlot[
          other1 /. Indeterminate -> 0]; 
ListPlot[
          other2 /. Indeterminate -> 0];
]

PlotLyapunov[0.5, 20, 200, 0.1, 30, 10]

PlotLyapunov[0.5, 10, 100, 0.1, 30, 720]

ListPlot[other /. Indeterminate -> 0];

ListPlot[other /. Indeterminate -> 0];