N1 <- c(1:200) N2 <- c(1:200) N3 <- c(1:200) yr <- c(1:200) cc <- c(0.01,0.01,0.01) a <- c(0.2,0.4,0.1,0.1,0.2,0.4,0.4,0.1,0.2) r <- c(0.5,0.6,0.7) rand1 <- runif(200,-0.5,0.5) rand2 <- runif(200,-0.5,0.5) rand3 <- runif(200,-0.5,0.5) for (i in 1:199) { N1[i+1] <- cc[1] + exp(r[1]+rand1[i]-a[1]*N1[i]-a[2]*N2[i]-a[3]*N3[i])*N1[i]; N2[i+1] <- cc[2] + exp(r[2]+rand2[i]-a[4]*N1[i]-a[5]*N2[i]-a[6]*N3[i])*N2[i]; N3[i+1] <- cc[3] + exp(r[3]+rand3[i]-a[7]*N1[i]-a[8]*N2[i]-a[9]*N3[i])*N3[i]; } op<-par(mfrow=c(3,1)) plot(yr,N1,col="red",type="o") plot(yr,N2,col="green",type="o") plot(yr,N3,col="cyan",type="o") par(op)