Fig_05

#---- Nicholson & Bailey

# variables
N <- numeric()
P <- numeric()

# coeficientes
lambda <- 2.8
c <- 1.5
a <- 0.009

# condiciones simulación
ngen <- 300      # generaciones
N[1] <- 150     # N inicial
P[1] <- 150     # P inicial

# simulación
for (t in 1:ngen) {
  f <- exp(-a * P[t])
  N[t + 1] <- lambda * N[t] * f
  P[t + 1] <- lambda * c * N[t] * a * P[t] * f
}

# equilibrio
Peq <- log(lambda) / a
Neq <- 1 / (a *c)

# gráfico
pdf("fig_05.pdf", width = 7.5, height = 6)
par(pty = "s",
    bty = "n",
    mar = c(5, 2, 1, 1),
    xaxs = "i",
    yaxs = "i",
    las = 1)
plot(N, P,
    type = "l",
    lty = 2,
    col = "lightblue",
    xlim = c(0, 200),
    ylim = c(0, 300))
points(N, P,
    cex = 0.8,
    pch = 16,
    col = "royalblue")
abline(h = Peq,
    col = "orange",
    lwd = 0.5)
abline(v = Neq,
    col = "orange",
    lwd = 0.5)
dev.off()