I’m pretty new to this and I’m trying to figure out how I can graph multiple results of a randomly generated number between (0:0.04) and I have it graphed but I’m so confused as to why the blue graph lines are angling back towards the initial start of the graph? I was expecting it to look more fan-like and not like it should be looping back towards the beginning?

library(deSolve)  

sir_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dB <- -M*B_cells - beta*Cb*B_cells - beta_2*Ct*B_cells + (g1*(Cb+Ct)/(g2+(Cb+Ct))) 
    dCb <- M*B_cells + beta*Cb*B_cells + beta_2*Ct*B_cells - alpha*Cb
    dT <- -nu_A*Cb*T_cells - nu_b*Ct*T_cells + (h1*(Cb+Ct)/(h2+(Cb+Ct)))
    dAt <- nu_A*Cb*T_cells + nu_b*Ct*T_cells - beta_2*Ct*At - beta*Cb*At
    dLt <- theta*(beta_2*Ct*At + beta*Cb*At)
    dCt <- (1-theta)*(beta_2*Ct*At + beta*Cb*At) - alpha_2*Ct
    dZ <- mu*Lt
    df <- -nu_F*Lt*f
    dIf <- nu_F*Lt*f
    return(list(c(dB, dCb, dT, dAt, dLt, dCt, dZ, df,dIf)))
  })
} 
mu_vec <- runif(5, min=0, max = 0.04)
res <- vector(length(mu_vec), mode = "list")
for (k in seq_along(mu_vec)) {
  res[[k]] <- ode(y = initial_values, times = time_values,
                  func = sir_equations, parms= c(mu=mu_vec[k],
                                                 M=0.05,
                                                 beta = 4.819e-4,
                                                 beta_2 =5e-4,
                                                 nu_A = 0.005,
                                                 nu_b = 0.01,
                                                 nu_F = 0.7,
                                                 alpha = 0.015,
                                                 alpha_2 = 0.010,
                                                 theta = 0.5,
                                                 g1 = 0.05,
                                                 g2 =0.001,
                                                 h1 = 0,
                                                 h2 = 10))
}  

names(res) <- mu_vec
dd <- dplyr::bind_rows(lapply(res, as.data.frame), .id= "mu")

dd$mu <- as.numeric(dd$mu)

initial_values <- c( 
  B_cells = 50  
  , Cb = 0 
  , T_cells=50 
  , At = 0
  , Lt = 0
  , Ct = 0
  , Z = 0
  , f =5 
  , If =0
  
)  


time_values <- seq(0,300) # hours


with(dd, {
  plot(x=time, y=B_cells, col="black", type="l", ylim= c(0,100),  xlab="Time (Hours)", ylab="Population density", main="Marek's Model")          #Plot the data for S over time
  lines(time, T_cells, col="red")  
  lines(time, Cb, col= "green")                                  #Add a line for I over time
  lines(time, At, col="blue")
  lines(time,Lt, col="purple" ) 
  lines(time,Ct, col="yellow" )  
  lines(time, Z, col= "lightblue") 
  lines(time, f, col="aquamarine") 
  lines(time,If, col="pink" ) 
})

I would greatly appreciate anyones input! Thanks so much!
Graph with weird lines

I have tried using the replicate function instead to see if maybe I should try another way but I got confused with how to implement the random values for mu and replicate function doesn’t seem to capture what I want to do..