library(deSolve)
library(ggplot2)
library(latex2exp)

sir<-function(t, init, parms) {
  with(as.list(c(init, parms)), {
    # nondimensional form
    ds <- -R0*s*i
    di <- R0*s*i - i
    list(c(ds, di))
  }) 
}

s0 <- 1 - 10^(-7)
T <- 100
init <- c(s=s0, i=1-s0)
times <- seq(0, T, by=.01)

df <- function(R0) {
  parms <- c(R0=R0)
  out <- ode(init, times, sir, parms)
  df <- data.frame(t=out[,1], i=out[,3], R0=R0)
  df
}
data <- do.call(rbind, lapply(seq(1.0, 3.0, by=.2), df))

p <- ggplot(data, aes(color=sprintf("%1.1f", R0))) +
  geom_path(aes(x=t, y=i)) + 
  scale_x_continuous(limits=c(0, T), breaks=seq(0, T, by=10)) +
  scale_y_continuous(limits=c(0, .31), breaks=seq(0, .31, by=.05)) +
  labs(
    x=TeX("$\\gamma\\mathit{t}$"),
    y=TeX("$\\mathit{I}/\\mathit{N}$"),
    color=TeX("$\\beta\\mathit{N}/\\gamma$")
  ) +
  scale_fill_brewer("Spectral") +
  theme_minimal(base_family = "Times")

svg(filename="infected.svg", width=16, height=9)
plot(p)
dev.off()