library(deSolve)

# Define your differential equation function
growth_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    # Calculate the time-varying growth rate (r) using an external function or any other source
    r <- get_time_varying_growth_rate(time)

    # Differential equation for population growth (logistic growth)
    dN_dt <- r * N * (1 - N / K)

    # Return the derivatives
    return(list(dN_dt))
  })
}

# Function to define the time-varying growth rate
get_time_varying_growth_rate <- function(time) {
  # In this example, we'll use a simple sinusoidal function as the time-varying growth rate
  # You can define any function based on your specific use case
  r_max <- 0.2  # Maximum growth rate
  r_min <- 0.05  # Minimum growth rate
  period <- 5  # Period of the sinusoidal function
  
  r <- (r_max + r_min) / 2 + (r_max - r_min) / 2 * sin(2 * pi * time / period)
  return(r)
}

# Define time points for the simulation
times <- seq(0, 20, by = 0.1)

# Initial state (initial population)
initial_state <- c(N = 10)

# Model parameters
parameters <- c(K = 100)  # Carrying capacity

# Solve the differential equation
output <- ode(y = initial_state, times = times, func = growth_model, parms = parameters)

# Plot the population growth over time
plot(output, xlab = "Time", ylab = "Population", main = "Population Growth Over Time")