Как я могу навязать начальные условия дискретному PDE при использовании deSolve?

Question

Я пытаюсь решить систему дифференциальных уравнений с одним ОДУ и одним УЧП второго порядка. Я дискретизировал PDE, чтобы решить его с помощью deSolve ode.1D, но теперь я не могу понять, как включить мои граничные условия для пространственной производной dE / dx. Существуют граничные значения для dE / dx при x = 0 и x = 1, которые важны для получения правильного решения. Как включить их в модель?

Мой код:

      BVmod1D <- function(time, state, parms, N, dx) {
  with(as.list(parms), {
    E <- state[1 : N]
    U <- state[(N + 1) : (2 * N)]
    

    coef1 <- Sv * io / (Qox - Qred)
    coef2 <- Tau * io / Cd
    BV <- exp(aa * f * (E - 0.5 * (1 + U) )) - exp(-ac * f * (E - 0.5 * (1 + U)))


    dEdx <- diff(c(E, E[N])) / dx # first spatial derivative of E
    ddEddx <- diff(c(dEdx, dEdx[N])) / dx # second spatial derivative of E
    

    dU <- coef1 * BV # dU/dt, Eqn 8
    dE <- (ddEddx - (coef2 * BV)) / Tau #  dE/dt, Eqn 12

    return(list(c(dE, dU)))
  })
}


pars <- c(Sv = 1.72e5, 
          io = 1.71e-6, 
          Qox = 1.56, 
          Qred = 0,
          Tau = 3.79e-6, 
          Cd = 7.42e-8, 
          aa = 0.7674, 
          ac = 0.2326, 
          ks = 0.042992, 
          sig = 1.67e-6, 
          f = 38.92237)


R <- 1
N <- 50
dx <- R / N
Vo <- 0.5


# Initial and Boundary Conditions

yini <- rep(0, (2 * N))
yini[ 1 : N ] <- 2 * Vo
yini[ (N + 1) : (2 * N)] 
times <- seq(0, 1 , by = 0.002 ) 


tail(out <- ode.1D(
  y = yini,
  times = times,
  func = BVmod1D,
  parms = pars,
  nspec = 2,
  N = N,
  dx = dx
))