R中的Desolve：返回的导数数错误

我目前正在try 创建一个函数，该函数将使用Desolve函数来解析舱模型.参数随着时间的推移而变化，基于温度和降雨量，所以我使用了ApxFun()来对参数估计进行内插.我的虚拟数据集中有4个变量，目前有122个时间步长.当我try 运行该函数时，我收到以下错误:

Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (488) must equal the length of the initial conditions vector (4)

我猜这与插补参数有关，因为很明显，4个变量*122个时间步长是488，但我真的很难弄清楚如何解决这个问题.

我将在下面粘贴我的代码！如有任何帮助，我们不胜感激.

#### Fake data set up ##########################################################

fake.dates <- seq(as.Date("2020/1/1"), as.Date("2020/5/1"), 'days')
# 122 dates 

fake.temp <- sample(-3:30, size = 122, replace = T)
all.precip <- seq(0, 30, by = 0.1)
fake.precip <- sample(all.precip, size = 122, replace = T)

fake.climate <- data.frame(fake.dates, fake.temp, fake.precip)

ost.sim.constanteggs <- function (E0, L0, L3f0, L3p0, start, end, temp, precip) {
  
  require(deSolve)
  
  ## Create a sequence of dates between start and end
  date.range <- seq(as.Date(start), as.Date(end), 'days')
  
  ## Get input time points (numeric)
  global.t <- seq(1, length(date.range))
  
  ## Parameters 
  dev = pmin(1, pmax(0, -0.07258 + (0.00976 * temp)))
  mig.1 = pmax(0, 0.0066 * precip - 0.0382)
  mig.2 = pmax(0, exp(-5.48240 + (0.45392 * temp) - (0.00540 * temp^2)))
  mu.1 = pmin(1, exp(-4.38278 - (0.10640 * temp) + (0.00540 * temp^2)))
  mu.2 = pmin(1, exp(-4.38278 - (0.10640 * temp) + (0.00540 * temp^2)))
  mu.3 = pmin(1, exp((-4.864 * temp) + (0.0048 * temp^2) + (0.00008 * temp^3)))
  mu.4 = pmin(1, exp((-5.743 * temp) + (0.0068 * temp^2) + (0.00003 * temp^3)))
  mu.5 = pmin(1, 10 * exp(-6.388 - (0.2681 * temp) + (0.01633 * temp^2) - (0.00016 * temp^3)))
  
  ## Interpolate the rates
  devrate <- approxfun(dev, method = 'linear', rule = 2)
  mig1rate <- approxfun(mig.1, method = 'linear', rule = 2)
  mig2rate <- approxfun(mig.2, method = 'linear', rule = 2)
  mu1rate <- approxfun(mu.1, method = 'linear', rule = 2)
  mu2rate <- approxfun(mu.2, method = 'linear', rule = 2)
  mu3rate <- approxfun(mu.3, method = 'linear', rule = 2)
  mu4rate <- approxfun(mu.4, method = 'linear', rule = 2)
  mu5rate <- approxfun(mu.5, method = 'linear', rule = 2)
  
  print("Parameters loaded successfully")
  
  # Model function 
  para.dyn <- function(times, para.init, para.par) {
    
    with(as.list(c(para.init, para.par)), {
      
      # Parameters
      dev1 <- devrate(times)
      mig1 <- mig1rate(times)
      mig2 <- mig2rate(times)
      mu1 <- mu1rate(times)
      mu2 <- mu2rate(times)
      mu3 <- mu3rate(times)
      mu4 <- mu4rate(times)
      mu5 <- mu5rate(times)
      
      # Differential equations 
      dE <- -(mu.1 + (2 * dev)) * E + 100
      dL <- -(mu.2 + (2 * dev)) * L + (2 * dev) * E
      dL3f <- -(mu.3 + mig.1) * L3f + (2 * dev) * L
      dL3p <- -mu.4 * (L3p * (1 - mig.2)) - mu.5 * (mig.2 * L3p) + mig.1 * L3f
      
      res <- c(dE, dL, dL3f, dL3p)
      list(res)
    })
  }
  
  print('Model function loaded successfully')
  
  para.init <- c(E = E0, L = L0, L3f = L3f0, L3p = L3p0)
  print('Inital conditions for state variables loaded successfully')
  
  para.output <- lsoda(y = para.init, 
                       times = global.t, 
                       func = para.dyn,
                       parms = NULL)
  
  return(para.output)
}


ss = "2020/1/1"
ee = "2020/5/1"

test <- ost.sim.constanteggs(E0 = 100, L0 = 0, L3f0 = 0, L3p0 = 0, start = ss, end = ee, temp = fake.climate$fake.temp, precip = fake.climate$fake.precip)

我希望得到的输出是每个变量(E、L、L3f、L3p)在每个时间点的值的数据帧.