Is there a different efficient way to link two models in R that use the deSolve library for calculations?

Question

I'm using the deSolve library in R to solve model_1. The output of model_1 is then used to update parameters in model_2. I run model_2 for one iteration, modify the parameters, and run it again. However, this loop takes a considerable amount of time due to the size of my model. Do you know of any alternative methods to speed up the simulation process? the r code is given below

library("deSolve")

#Define the model function`
dq <- function(t, x, parms) {
  dim <- parms[1]  
  parms <- parms[-1] 
  compsys <- matrix(parms[1:(dim * dim)], nrow = dim)  
  qout <- x %*% compsys  # calculate dq/dt
  list(qout)  
}
solve_eqn <- function(transfers, compnames, initial, times) {
  dim <- length(transfers[1,])  
  k <- transfers 
  idx <- c(1:dim)
  for (i in idx) {
    k[i, i] <- (-sum(transfers[i, ]))
  }
  odesolution <- radau(func = dq, y = initial, times = times, parms = c(dim, k))
  solution <- matrix(data = odesolution[, -1], nrow = length(times), byrow = FALSE)
  colnames(solution) <- compnames
  return(cbind(times,solution))
}

# define model_1
model_1 <- matrix(c(0, 2, 0.2,
                 0.1, 0, 1,
                 0.8, 0, 0),
               nrow = 3, ncol = 3, byrow = TRUE)

colnames(model_1) <- c("blood", "ECFs", "ECFf")
rownames(model_1) <- c("blood", "ECFs", "ECFf")

#define model_1 input
initial <- c(blood = 1000, ECFs = 0, ECFf = 0)
interval <- seq(1, 10, 1)

#solve model_1
result <- solve_eqn(transfers = model_1, compnames = colnames(model_1), initial = initial, times = interval)

# define model_2
model_2 <- matrix(c(1, 2, 2,
                 4, 0, 1,
                 0.8, 0, 3),
               nrow = 3, ncol = 3, byrow = TRUE)
colnames(model_2) <- c("kidney", "liver", "stomach")
rownames(model_2) <- c("kidney", "liver", "stomach")

#define model_2 input
initial2 <- c(kidney=1000, liver=0, stomach=0)

decay_in_comps <- NULL
for (t in interval)
{
  #check if the quantity at time t exists in result
  #if yes, replace the element of data2 with quantity 
  row_number <- which(result["times"]==t)
  if (length(row_number) != 0)
  {
    model_2["kidney", "liver"] <- as.numeric(result[row_number,"liver"])
    model_2["stomach","kidney"] <- as.numeric(result[row_number,"kidney"])

  }
  #solve the model_2  for one iteration using radau function
  result2 <- solve_eqn(transfers = model_2, compnames = colnames(model_2), initial = initial2, times = interval)
  #concatenate the results of each iteration
  decay_in_comps <-rbind(decay_in_comps,initial2)
  rownames(decay_in_comps) <- NULL
  #change the input parameters for solver with each iteration
  initial2 <-result2[2,-1]
}
print(cbind(interval,decay_in_comps))