CL <- 5
V1 <- 19
Rate <- 1000
Q <- 0.05
V2 <- 5.263157895
time <- seq(0, 48, 1)

library(deSolve)
library(dplyr)

#deSolve function
comp2Inf <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    #Differential equations
    dC1 <- Rate/V1 - C1*(CL/V1) - C1*(Q/V1) + C2*(Q/V1)
    dC2 <- C1*(Q/V2) - C2*(Q/V2)
    #return rate of change
    list(c(dC1, dC2))
  })
}

ref_fcn <- function(times, Rate, CL, V1, Q, V2) {
  b = (Q/V1+Q/V2+CL/V1-sqrt((Q/V1+Q/V2+CL/V1)**2-4*Q*CL/V1/V2))/2
  a = Q*CL/V1/V2/b
  A = (a-Q/V2)/V1/(a-b)
  B = (b-Q/V2)/V1/(b-a)
  return(sapply(times, function(t) {
    Rate*(A/a*(1-exp(-a*t))+B/b*(1-exp(-b*t)))
  }))
}

my_fcn <- function(times, Rate, CL, V1, Q, V2) {
  a = (Q/V2+Q/V1+CL/V1-sqrt((Q/V2+Q/V1+CL/V1)**2-4*CL*Q/V1/V2))/2
  b = (Q/V2+Q/V1+CL/V1+sqrt((Q/V2+Q/V1+CL/V1)**2-4*CL*Q/V1/V2))/2
  return(sapply(times, function(t) {
    Rate*(Q/a/b+(Q-(V2*a))*exp(-a*t)/a/(a-b)+(Q-(V2*b))*exp(-b*t)/b/(b-a))
  }))
}

tab_deS <- ode(y = c(C1=0,C2=0), times = time, func = comp2Inf, parms = list(Rate=Rate, CL=CL, V1=V1, Q=Q, V2=V2)) %>% as.matrix
arr_ref <- ref_fcn(time, Rate, CL, V1, Q, V2)
arr_my <- my_fcn(time, Rate, CL, V1, Q, V2)

data.frame(time=time, deS=tab_deS[,2], ref=arr_ref, my=arr_my)