may = function(startr,endr,stepr, k, noise, n0, nostep,singleplot,animation) {
# animated plotting of a time series of the May's equation (if singleplot>0)
# plus the final plot of N against r (bifurcation plot)

#startr,endr,stepr - minimum, maximum growth rate and difference between adjacent values
#k - carrying capacity
#n0 - initial population size
#nostep - number of steps for simulation

# singleplot == 0: no plot except the bifurcation diagram
# singleplot == 1: N against time
# singleplot == 2: state space
# animation - whether to animate

  startstep =min(50,nostep-1)
  if (startr >= endr | stepr<=0) return("Error");

  dots=NULL

  for (r in seq(startr,endr, by=stepr)) {
     a=maysingle(r,k,noise,n0,nostep);
     dots=cbind(dots, a[startstep:nostep])
     if (singleplot==1) {
        plot(1:(nostep+1),a, xlim=c(1,nostep+1), ylim=c(0, k*3/2), type="l",lwd=2,col="blue",main=paste("R value:",r),
             xlab="Step", ylab="Population size");
        Sys.sleep(animation);
     }
     if (singleplot==2) {
        parabolasingle(r,k,noise,n0,nostep,endr);
        Sys.sleep(animation);
     }
   }
   dev.new(3)
   matplot(seq(startr,endr, by=stepr),t(dots),type="p",pch=20,col="black",xlab="R value",ylab="Final population sizes")
}