Speed up parallelised solving of ODE containing events

Question

The following code represents a simplified version of my actual problem:

# Libraries
library(desolve)
library(parallel)

# End time
max_time <- 300

# Time vector
time <- seq(0, max_time, 0.1)

# Number of events
n_events <- c(6, 60, 600, 6000)

# p1 values
p1_vars <- 10^c(-4, -2, 0, 2)

# Parameters
p <- data.frame(p1 = 1.14*24,
                p2 = 0.14*24,
                p3 = 0.06*24,
                p4 = 0.5,
                p5 = 0.05,
                p6 = 0,
                p7 = 0.218,
                p8 = 0.477,
                p9 = 0.5/6)

# Initial conditions
x0 <- c(x1 = 0,
        x2 = 0,
        x3 = 2e11,
        x4 = 8e11)

# Model
sys <- function(t, y, p) {
  dy = numeric(4)

  dy[1] = p$p2*y[2] - (p$p3 + p$p1)*y[1]
  dy[2] = p$p3*y[1] - p$p2*y[2]
  dy[3] = (p$p4 - p$p7 - p$p8)*y[3] + p$p5*y[4] - p$p9*y[1]*y[3]
  dy[4] = p$p7*y[3] - (p$p5 + p$p6)*y[4]
  list(dy)
}

# Number of corse to parallelise
n_cores <- detectCores() - 1

# Initialise
input.df <- data.frame(eventfreq = rep(n_events/max_time, length(p1_vars)),
                       par = rep(p1_vars, each = length(n_events)))

res.df <- data.frame(eventfreq = NaN,
                     par = NaN,
                     out = NA)

# Parallel function
par_fcn <- function(i) {

  # p1 value
  p$p1 <- input.df$par[i]

  # Event value
  ev_val <- 1000/input.df$eventfreq[i]*max_time

  # Events
  ev <- data.frame(var = 'x1',
                   time = seq(0, max_time, by = 1/input.df$eventfreq[i]),
                   value = 100,
                   method = 'add')

  # Solve ODE
  sim.df <- as.data.frame(ode(y = x0, times = time, func = sys, p = p, events = list(data = ev)))

  # Fill
  res.df$eventfreq <- input.df$eventfreq[i]
  res.df$par <- input.df$par[i]
  res.df$out <- sim.df$x3 + sim.df$x4

  # Return
  return(res.df)
}

# Evaluate in parallel
res_parallel <- mclapply(seq_len(nrow(input.df)), par_fcn, mc.cores = n_cores)
out.df <- do.call(rbind, res_parallel)

As you can see, I try to solve an ODE with desolve for a range of parameter values and different number of events. Even though, I parallelised the code and spread it across 7 cores, my actual problem still takes very long. I've identified the ODE solver as the bottleneck, especially when the number of events gets large.

Are there any ideas how this code could be speeded up?

I've also been looking into implementing this in compiled code according to the vignette R Package deSolve: Writing Code in Compiled Languages. But I'm too inexperienced in C to do it completely on my own.

Here's what I got so far:

/*file model.c */
#include <R.h>
static double p[9];
#define p1 p[0]
#define p2 p[1]
#define p3 p[2]
#define p4 p[3]
#define p5 p[4]
#define p6 p[5]
#define p7 p[6]
#define p8 p[7]
#define p9 p[8]

/* initialiser */
void initmod(void (*odeparms)(int *, double *))
{
    int N = 9;
    odeparms(&N, p);
}

/* Events */
void event(int *n, double *t, double *y)
{
    y[0] = y[0] + 100;
}

/* Derivatives */
void derivs (int *neq, double *t, double *y, double *ydot, double *yout, int *ip)
{
    if (ip[0] < 1) error("nout should be at least 1");
    ydot[0] = p2*y[1] - (p3 + p1)*y[0];
    ydot[1] = p3*y[0] - p2*y[1];
    ydot[2] = (p4 - p7 - p8)*y[2] + p5*y[3] - p9*y[0]*y[2];
    ydot[3] = p7*y[2] - (p5 + p6)*y[3];

    yout[0] = y[0];
    yout[1] = y[2] + y[3];
}

/* END file model.c */

Therein, I especially can't follow the vignette of how to pass the n_events and p1_vars parameter to the C code.

Any help or hints would be greatly appreciated.