How to force boundary conditions on discretized derivative in deSolve?

Hello,

I am using R package deSolve to solve a system of two differential
equations for a one-dimensional spatial and time-based problem. There is
one ODE and a second-order PDE. In order to solve with the function ode.1D,
I've discretized the spatial derivative and put both equations in terms of
the time derivative only. However, I'm now stuck with the problem of being
unable to impose boundary conditions on the spatial derivative for flux at
the edges of the system. How can I force a value for the discretized dE/dx
part of my equations at x = 0 and x = 1?

The code I have been using is below. The ode.1D solver will run on it, but
the solutions aren't correct for my system owing to the missing boundary
conditions.

Thanks,

Molly C Brockway
Materials Science - PhD Candidate
Metallurgical and Materials Engineering - B.S.

Montana Technological University


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

    coef1 <- Sv * io / (Qox - Qred)
    coef2 <- Tau * io / Cd
    BV <- exp(aa * f * (E - 0.5 * (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
    dE <- (ddEddx - (coef2 * BV)) / Tau

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

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)] <- 1
times <- seq(0, 1 , by = 0.002 )

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