I am continuing in Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Lifeby Steven Strogatz. To get a feeling on it, I was building a group of things which have only a minute influence on each other are able to synchronize their behavior. I got some more explanations how things fit together and adapted my calculations. Two things are changed.
You can actually use one of the items as a kind of clock, only register status when it goes to zero. This makes it possible to plot behavior over much longer times. Then, I found my function of last week a bit primitive, by going step-wise through time. I was thinking at first to reduce the step-size  why not do the more exact thing and follow a function? So, I scaled the logit function a bit to provide behavior of an item and build some functions to determine when an item goes over the threshold. Code is at the bottom.

Even with a small additive effect they synchronize, as shown below. They are not completely synchronized, but close enough for me.

R code

library(ggplot2)
library(arm)


plotpart <- function(xmat,fname='een.png',thin=1) {
  nstep <- ncol(xmat)
  coluse <- seq(1,nstep,by=thin)
  xmat <- xmat[,coluse]
  fin <- length(coluse)
  nitems <- nrow(xmat)
  df <- data.frame(score=as.numeric(xmat),
      cycle=rep((coluse),each=nitems),
      item =rep(1:nitems,fin))
  g<- ggplot(df,aes(x=cycle,y=score,group=item,alpha=score)) + 
      geom_line() + 
      scale_alpha_continuous(range=c(.02,.0201) )+
      theme(legend.position=’none’) +
      xlab(‘Iteration’) + theme(panel.background = element_rect(fill=’white’))
  png(fname)
  print(g)
  dev.off()
}
time2score <- function(x) 2*invlogit(x)-1
score2time <- function(x) logit((1+x)/2)

time_delta_score <- function(score1,score2=.99) {
  score2time(score2)-score2time(score1)
}
score_delta_time <- function(scorenow,tdelta) {
  tnow <- score2time(scorenow)
  time2score(tnow+tdelta)
}

onestep <- function(score,limit=0.99,spil=1e-4) {
  mscore <- max(score)
  dtnew <- time_delta_score(mscore,limit)
  scorenew <- score_delta_time(score,dtnew)
  maxed1 <- maxed <- scorenew==max(scorenew)#>limit*.99999999
  while(sum(maxed)>0) {
    scorenew[!maxed] <- scorenew[!maxed] + sum(maxed)*spil
    scorenew[maxed] <- 0
    maxed <- scorenew>limit
    maxed1 <- maxed | maxed1
  }
  scorenew[maxed1] <- 0
  scorenew
}

nitems <- 500
niter <- 10000
xmat <- matrix(0,nrow=nitems,ncol=niter)

score  <- time2score(runif(nitems,0,score2time(.99)))
n <- 0
while (n < niter) {
  score <- onestep(score,spil=1e-7)
  if (score[1]==0) {
    n<- n+1
    xmat[,n] <- score
  }
}

plotpart(xmat[1:250,],’een.png’,thin=20)