dt = 0.0001                           # time step, delta t
tmax = 25                             # max time, s
t = np.linspace(0, tmax, num=tmax/dt) # time vector
nt = len(t)                           # total number of time steps

def rates(c, t):
    """
    w = wood-oil as conc[0]
    nv = non-volatiles as conc[1]
    v = volatiles as conc[2]
    """
    Knv = 0.3
    Kv = 0.8
    rw = -(Knv + Kv) * c[0]
    rnv = Knv * c[0]
    rv = Kv * c[0]
    return [rw, rnv, rv]

cc = sp.odeint(rates, [1, 0, 0], t)

def dCdt(c, t):
    """
    w = wood-oil as conc[0]
    nv = non-volatiles as conc[1]
    v = volatiles as conc[2]
    """
    Knv = 0.3
    Kv = 0.8
    Kr1 = 9.2e7;    r1 = 2.5
    Kcr = 4.1e-5;   cr = 0.9
    Kc1 = 3.7e5;    c1 = 1.1
    Kd = 8.0e-4;    d = 0.9
    Ka = 6.1e-6;    a = 1.4
    Kg = 1.8e-4;    g = 0.8
    Kr2 = 37.0e5;   r2 = 0.7
    rw = -(Knv + Kv)*c[0]
    rnv = Knv*c[0] - Kr1*c[1]**r1 - Kcr*c[1]**cr - Kc1*c[1]**c1
    rv = Kv*c[0] + Kcr*c[1]**cr - Kd*c[2]**d - Ka*c[2]**a - Kg*c[2]**g - Kr2*c[2]**r2
    return [rw, rnv, rv]

cc2 = sp.odeint(dCdt, [1, 0, 0], t)