# Load the deSolve package
library(deSolve)

# Define the system dynamics model
model <- function(time, state, parameters) {
  
  # Extract the variables from the state vector
  S <- state[1] # Susceptible population
  I <- state[2] # Infected population
  R <- state[3] # Recovered population
  
  # Extract the parameters
  beta <- parameters[1] # Infection rate
  gamma <- parameters[2] # Recovery rate
  
  # Define the differential equations
  dS <- -beta * S * I          # Susceptible population decreases
  dI <- beta * S * I - gamma * I  # Infected population increases then decreases
  dR <- gamma * I             # Recovered population increases
  
  # Return the rates of change
  return(list(c(dS, dI, dR)))
}

# Define the initial conditions
state <- c(S = 999, I = 1, R = 0)

# Define the parameters
parameters <- c(beta = 0.3, gamma = 0.1)

# Define the time steps
times <- seq(0, 100, by = 1)

# Solve the system dynamics model
output <- ode(y = state, times = times, func = model, parms = parameters)

# Plot the output
plot(output, xlab = "Time", ylab = "Population", main = "Spread of Virus")
legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = 1:3, lty = 1)