Solving ODEs in R language, as example for test this class:

 func <- function(t, x, parms) {
    with(as.list(c(parms, x)), {
    
       dP <- a * P      - b * C * P
       dC <- b * P * C  - c * C
      
       list(c(dP, dC))
    })
 }

 y0    <- c(P = 2, C = 1)
 parms <- c(a = 0.1, b = 0.1, c = 0.1)
 out   <- ode(y = y0, times = 0:100, func, parms = parms)
 
 head(out)
 plot(out)