这里是一个解决方案的一部分。时间被四舍五入与floor和precip转换为一个强制函数与常数插值。
igrow <- ifelse(...)-函数看起来很奇怪，尤其是模型中的值永远不应该是NA。

library(deSolve)

mod <- function(t, yini, param, f_precip) {   
  app1 <- ifelse(floor(t) == 156, 1.12, 0)   # use round, floor or trunc
  jgrow <- 32
  
  ## the following can also be converted to a forcing function
  ## but I suspect it may have some errors
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t > 168, 0.9, ifelse(t < 136, 0, NA)))
  
  cover <- 0.9 * igrow / jgrow
  P <- f_precip(t)
  
  with(as.list(c(yini, param)), { 
    mharv <- ifelse(floor(t) == 201, mf, 0) # use floor, trunc or round
    d_MF <- (app1 * sa * (1 - drift) * cover * sa) - (mf * exp(-kf * dt)) + 
            (yf * (mf * exp(-kf * dt))) - mf * (1 - exp(-fet * P * dt)) - mharv
    return(list(c(d_MF))) 
  })
} 

params <- list(drift = 0.1, dt = 1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)

t <- seq(1, 365, by = 1)
Precip <- rep(0, 365)
Precip[155:205] <- c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 
                     0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 
                     0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 
                     3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)

## define a forcing function by tabular interpolation
f_precip <- approxfun(t, Precip, method = "constant", rule = 2)

y0 <- c(mf = 0)
sol <- ode(y = y0, times = t, func = mod, parms = params, 
           f_precip = f_precip)

plot(sol)

编辑

为了回答OP的评论，为什么“应用后第2天功能达到最大值”，我们可以应用以下修改：
在模型函数中，将app1替换如下。在状态向量c()之后添加一些内部变量作为return中的额外输出也是一个好主意。

mod <- function(t, yini, param, f_precip, f_app) {   
  #app1 <- ifelse(floor(t) == 156, 1.12, 0)
  app1 <- f_app(t)
    .....
    return(list(c(d_MF), P = P, app = app1, igrow=igrow))
  })
}

然后定义一个新的强制：

f_app <- approxfun(c(0, 156,        156+1/24, 365),
                   c(0, 1.12 * 24,    0,        0), method="constant", rule=2)

我们看到，只需要相关的时间步长。
然后必须修改对ode的调用以包含新的强制。减少时间步长非常重要，这样求解器的自动步长算法就不会忽略外部强制。可以通过在t-向量中设置较小的时间步长或通过附加参数hmax来完成。

sol <- ode(y = y0, times = t, func = mod, parms = params, 
           f_precip = f_precip, f_app = f_app, hmax=1/24)

plot(sol, mfrow=c(2, 2))

然而，一个小的时间步长会减慢模拟的速度。为了避免这种情况，可以使用事件机制。然后，您必须重新制定模型，使其在导数中不包含app1，而是计算它对mf的影响。然后，在事件发生时，它被直接和瞬间添加到mf。

app1和mharv实际上始终为零，几乎可以肯定斜率函数不会在该特定点处进行评估，即使如此，该单点尖峰也仅具有可忽略的影响。
出于同样的原因，在调用mod时，t几乎肯定不会是整数，当试图使用t作为整数索引时，应该会得到一个错误，但可能会自动截断。
所以最后，导数表达式中的每一项都有一个因子为零，使零函数成为正确的解。