library(deSolve); #Parameters parms <- c(F_IIF = 0.1, k_IFP = 0.05, F_IDF = 0.3, k_DFSC = 0.05, k_DFP = 0.05, F_ISC = 0.05, k_SCD = 0.5, k_SCP = 0.2) #Differential equations Rumen <- function(t, q, parms) { with(as.list(c(parms, q)), { dQ_IF <- F_IIF - k_IFP*Q_IF dQ_DF <- F_IDF - k_DFSC*Q_DF - k_DFP*Q_DF dQ_SC <- F_ISC + k_DFSC*Q_DF - k_SCD*Q_SC - k_SCP*Q_SC res <- c(dQ_IF, dQ_DF, dQ_SC) list(res) }) } #Initial conditions qinit <- c(Q_IF= 1, Q_DF = 1.5, Q_SC = 0.5) #Time steps times <- seq(0, 500, by = 0.2) #Numerical solution out1 <- rk(qinit, times, Rumen, parms, method = "rk4") #Figure with numerical solution plot(out1)