What is the problem with my R code for solving and plotting a coupled system of ODE's?

I'm reasonably new to R and have this system of ODEs
\frac{dS}{dT} = -\beta(1-\mu)S(t)\frac{I(t))}{1-D(t))}
\frac{dI}{dT} = -\beta(1-\mu)S(t)\frac{I(t))}{1-D(t))}-\delta I(t))-\phi I(t))
\frac{dR}{dT}=\delta I(t))
\frac{dD}{dT}=\phi I(t))

I'm trying to create a plot of S, I, R and D with respect to time and have written this code with my basic knowledge of R and the help of similar examples in the deSolve manual. However the plots my code generates are not what I was expecting.

(Note. the parameter values I am using are arbitrary, I'm planning to use a sums of squares approach to adjust these to the true parameter values in my dataset which is COVID-19 case and death data from Mexico, also I misnamed the parameter letters in my code, they represent the same values)

This is my code: 


```
install.packages("deSolve")
library("deSolve")

sird_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
dS <- -beta * (1- mu) * S * (I/(1-D))
dI <- -beta * (1- mu) * S * (I/(1-D))-sigma * I - theta * I
dR <- sigma * I 
dD <- theta * I 
return(list(c(dS, dI, dR, dD)))
})
} 

parameter_values <- c( 
beta = 0.95, #infectious contact rate (/person/day)
sigma = 0.5, #recovery rate (/day)
theta = 0.043, #mortality rate (/day)
mu = 0.1 #adjustment for efficacy of countermeasures
)

inital_values <- c(
S = 126000000, #suseptible population (adjust based on ylur data)
I = 1, #infected on day 1 of data 
R = 0, #recovered on day 1
D = 0 #dead on day 1
)


time_values <- seq(0, 10) #days(adjust based on your data)
ls()

sird_values_1 <- ode(
y = inital_values,
times = time_values,
func = sird_equations,
parms = parameter_values 
)

sird_values_1 <- as.data.frame(sird_values_1)
with(sird_values_1, {
plot(time, S, type = "l", col = "blue",
xlab = "time (days)", ylab = "# of people")
lines(time, I, col = "black")
lines(time, R, col = "red")
lines(time, D, col = "orange")
})

legend("right", c("susceptible", "infectious", "recovered", "dead"),
col = c("blue", "black", "red", "orange"), lty = 1, bty = "n")
```

I was expecting a plot like this But what I get is this

I have tried adjusting the parameter values but nothing I do has much impact on the plot and I would really appreciate support in understanding where I have gone wrong.

Here is a model based on the SIRD model described on wikipedia, with your mu adjustment:

sird_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- - (beta * (1 - mu) * S * I ) / (S + I + R + D)
    dI <- (beta * (1 - mu) * S * I) / (S + I + R + D) - sigma * I - theta * I
    dR <- sigma * I 
    dD <- theta * I 
    return(list(c(dS, dI, dR, dD)))
  })
} 

I changed the time steps:

time_values <- seq(0, 100, 0.01)

Then I also changed the plotting code, to take the correct y-limits:

with(sird_values_1, {
  plot(
    time, S, 
    ylim = c(0, max(sird_values_1[-1])),
    type = "l", col = "blue",  xlab = "time (days)", ylab = "# of people"
  )
  lines(time, I, col = "black")
  lines(time, R, col = "red")
  lines(time, D, col = "orange")
})

Or using ggplot2:

sird_values_1 |> 
  tidyr::pivot_longer(-time) |> 
  ggplot(aes(time, value, color = name)) +
  geom_line() +
  scale_color_manual(values = c(S = "blue", I = "black", R = "red", D = "orange")) +
  scale_y_continuous(labels = scales::label_number()) +
  labs(y = '# of people', x = 'time (days)', color = NULL) +
  theme_classic()