Func2, times, y, rho error

I keep getting the error:

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

I am trying to create a model based off of Brigatti et al 2009 (pred-prey model w a spatial component)

x<-c(1:40000)

left_shift = function(x) {
  x[c(2:length(x), 1)]
}

right_shift = function(x) {
  x[c(length(x), 1:(length(x) - 1))]
}

laplace = function(x) {
  return(c(left_shift(x) + right_shift(x) - 2 * x))
}

dxdt <- function(time, state, pars) {
  prey = state[1:length(state) / 2]
  pred = state[(length(state) / 2 + 1):length(state)]
  dprey = pars[5] * laplace(prey) + pars[1] * prey - x[2] * prey * pred
  dpred = pars[5] * laplace(pred) + pars[3] * prey * pred - pars[4] *  pred
  list(c(prey, pred, dprey, dpred))
}

time <- seq(0, 600, by = 1)
pars <- c(alpha=1, 
      beta = .5, 
      gamma = .2, 
      delta = .6,
      D = 0.000008 #(0.004*0.004/2), #diffusion coefficient
        )

state <- rep(0.1, 51)

out <- as.data.frame(ode(func = dxdt, y = state, parms = pars, times =  time))

Solution

A few problems. First, missing parentheses.

prey = state[1:length(state) / 2]

should read

prey = state[1:(length(state) / 2)]

Second, your initial conditions are an odd number in length. state should specify the initial conditions for both prey and predator (in that order). So, for each location there should be two values and, consequently, the vector should always be a multiple of two in length.

Thirdly, your function dxdt should return list(c(dprey, dpred)). There is no reason to return the values for the state variables, because the ODE solver will calculate those.

Fix those and this is what you get:

left_shift = function(x) {
  x[c(2:length(x), 1)]
}

right_shift = function(x) {
  x[c(length(x), 1:(length(x) - 1))]
}

laplace = function(x) {
  return(c(left_shift(x) + right_shift(x) - 2 * x))
}

dxdt <- function(time, state, pars) {
  prey = state[1:(length(state) / 2)]
  pred = state[(length(state) / 2 + 1):length(state)]
  dprey = pars[5] * laplace(prey) + pars[1] * prey - x[2] * prey * pred
  dpred = pars[5] * laplace(pred) + pars[3] * prey * pred - pars[4] *  pred

  list(c(dprey, dpred))
}

time <- seq(0, 600, by = 1)
pars <- c(alpha=1, 
          beta = .5, 
          gamma = .2, 
          delta = .6,
          D = 0.000008 #(0.004*0.004/2), #diffusion coefficient
)

state <- rep(0.1, 50)

out <- as.data.frame(ode(func = dxdt, y = state, parms = pars, times =  time))