library(deSolve)


model <- function(t, Y, parameters) {
  with (as.list(parameters),{

    dy1 = -k1*Y[1] + k2*Y[2]*Y[3]
    dy3 = k3*Y[2]*Y[2]
    dy2 = -dy1 - dy3

    list(c(dy1, dy2, dy3))
  })
}


jac <- function (t, Y, parameters) {
  with (as.list(parameters),{

    PD[1,1] <- -k1
    PD[1,2] <- k2*Y[3]
    PD[1,3] <- k2*Y[2]
    PD[2,1] <- k1
    PD[2,3] <- -PD[1,3]
    PD[3,2] <- k3*Y[2]
    PD[2,2] <- -PD[1,2] - PD[3,2]

    return(PD)
  })
}


parms <- c(k1 = 0.04, k2 = 1e4, k3 = 3e7)
Y     <- c(1.0, 0.0, 0.0)
times <- c(0, 0.4*10^(0:11))
PD    <- matrix(nrow = 3, ncol = 3, data = 0)

out   <- ode(Y, times, model, parms = parms, jacfunc = jac)

plot(out)


system("R CMD SHLIB mymod.c")

dyn.load(paste("mymod", .Platform$dynlib.ext, sep = ""))

## benchmarks

system.time(
  for(i in 1:10) out   <- ode(Y, times, model, parms = parms)
)

system.time(
  for(i in 1:10) out   <- ode(Y, times, model, parms = parms, jacfunc = jac)
)

system.time(
  for(i in 1:10)
    out <- ode(Y, times, func = "derivs", parms = parms,
               dllname = "mymod",
               initfunc = "initmod", nout = 1, outnames = "Sum")
)

system.time(
  for(i in 1:10)
    out <- ode(Y, times, func = "derivs", parms = parms,
               jacfunc = "jac", dllname = "mymod",
               initfunc = "initmod", nout = 1, outnames = "Sum")
)