Chunk con todos los elementos para modelar el caso.

library(deSolve)
InitialConditions <- c(crimenes.urb = 100,
                       crimenes.suburb = 50,
                       crimenes.rural = 200
  )

times <- seq(0,#initial time
             50,#end time
             1/4) #time step

politica <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    #Endogenous auxiliary variables
    crimenes.totales = crimenes.urb + crimenes.suburb + crimenes.rural
      
    #Flow variables
    incremento.crimenes.urb = poblacion.urb*tasa.incremento
    incremento.crimenes.suburb = poblacion.suburb*tasa.incremento
    incremento.crimenes.rural = poblacion.rural*tasa.incremento
    reduccion.crimenes.urb = ifelse(t>=años.politica.efectiva, tasa.reduccion.urb*crimenes.urb, 0)
    reduccion.crimenes.suburb = ifelse(t>=años.politica.efectiva, tasa.reduccion.suburb*crimenes.suburb, 0)
    reduccion.crimenes.rural = ifelse(t>=años.politica.efectiva, tasa.reduccion.rural*crimenes.rural, 0)
      
    #State (stock) variables (d)
    d.crimenes.urb = incremento.crimenes.urb - reduccion.crimenes.urb
    d.crimenes.suburb = incremento.crimenes.suburb - reduccion.crimenes.suburb
    d.crimenes.rural = incremento.crimenes.rural - reduccion.crimenes.rural
    
    list(c(d.crimenes.urb, d.crimenes.suburb, d.crimenes.rural), crimenes.totales = crimenes.totales)
  })
}

parameters<-c(poblacion.urb = 5000000,
              poblacion.suburb = 7000000,
              poblacion.rural = 4000000,
              tasa.incremento = 10/1000000,
              tasa.reduccion.urb = 5/100,
              tasa.reduccion.suburb = 7/100,
              tasa.reduccion.rural = 10/100,
              años.politica.efectiva = 10
  )


intg.method<-c("rk4")

out <- ode(y = InitialConditions,
           times = times,
           func = politica,
           parms = parameters,
           method =intg.method )

plot(out,
     col=c("blue"))

#Mejora de la salida de los gráficos.

library(ggplot2)
library(tidyr)


out_df <- as.data.frame(out)

# 1. Gráfico de crímenes por zona
out_long <- pivot_longer(out_df, 
                         cols = c(crimenes.urb, crimenes.suburb, crimenes.rural), 
                         names_to = "Zona", 
                         values_to = "Crimenes")

ggplot(out_long, aes(x = time, y = Crimenes, color = Zona)) +
  geom_line() +
  labs(title = "Crímenes por zona", x = "Tiempo", y = "Crímenes")

# 2. Gráfico de crímenes totales
ggplot(out_df, aes(x = time, y = crimenes.totales)) +
  geom_line(color = "red") +
  labs(title = "Crímenes totales", x = "Tiempo", y = "Crímenes totales")