I'm trying to derive the rate law of a first order followed by second order consecutive reaction:

$$ \frac{\mathrm{d}[\ce{A}]}{\mathrm{d}t}=-k_1[\ce{A}]$$ $$ \frac{\mathrm{d}[\ce{B}]}{\mathrm{d}t}=k_1[\ce{A}]-k_2{[\ce{B}]^2}$$ $$ [\ce{A}](0)=[\ce{A}]_0,\ce{[B]}(0)=0$$

The solution for A is trivial: $$ \ce{[A]}=\ce{[A]}_0\exp(-k_1t)$$

The problem is B, whose rate law ODE is not linear: $$ \frac{\mathrm{d}[\ce{B}]}{\mathrm{d}t}=k_1[\ce{A}]_0\exp(-k_1t)-k_2{[\ce{B}]^2}$$

This is a nonlinear first order ODE, but can be rewritten to a linear second order ODE by substituting $[\ce{B}]$ with: $$[\ce{B}]=\frac{\mathrm{d}y}{\mathrm{d}t}\frac{1}{k_2y}$$

Which gives: $$ -\frac{\mathrm{d}^2y}{\mathrm{d}t^2}+k_1k_2[\ce{A}]_0\exp(-k_1t)y=0$$

The solution is quite complicate and contains Bessel functions of both the first and the second kind:

$$ [B] = \frac{C_1\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,J_1\left(\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-\ce{A}_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\,\sqrt{-A_{0}\,k_{1}\,k_{2}}-\ce{C}_2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,Y_1\left(-\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{2}\,\left(C_1\,J_0\left(\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)+C_2\,Y_0\left(-\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\right)}$$

There are two arbitrary constants but I only have one condition: $$\ce{[B]}(0)=0$$

So my questions are:

  • Is my approach solving the rate law correct ?
  • If so, how do I get rid of both constants ?
  • Is there any approximated solution to this rate law ?

[update]: The ODE is solved numerically with Runge-Kutta method. The symbolic solution seems correct.

  • $\begingroup$ Unless you really want an algebraic solution you can numerically solve $dB/dt$ quite easily and you know both initial conditions; at $t=0, [A]=[A_0],[B]=0$. Once you make a 2nd order eqn, you will have 2 initial conditions making 3 in all. I'm not sure what this could be, $[B]=0, t=\infty$ is probably not helpful, so this means knowing $dB/dt$ at some time. Your substitution is a cunning one :), its a pity it does not make a Bernoulli eqn. $\endgroup$
    – porphyrin
    Jun 25, 2018 at 21:57
  • 1
    $\begingroup$ @porphyrin I don't see how numerical solution is easy because of the two parameters k1 and k2. To solve that one needs to nest the numerical ode solver inside optimization iterations. The optimization of k1 and k2 will be a global fit of a 2D array. $\endgroup$
    – 7E10FC9A
    Jun 26, 2018 at 15:26
  • 1
    $\begingroup$ Of course you have to give k1 and k2 values which you will also have to do for your algebraic solution when you want to plot it. The form of the chemical kinetics is that A will decay exponentially while B will rise and fall until no B is left but only the product P. This forms only from B. The product can be calculated by mass balance since $\ce{A0 + B0 =A + B + P}$ at all times. If you have real data to fit then it is still routine to do nonlinear optimisation with numerical solution of equations. $\endgroup$
    – porphyrin
    Jun 26, 2018 at 16:02

1 Answer 1


The following code is in Python and uses numpy, scipy and matplotlib that are usually packaged with python. The code will solve the equations $\mathrm{A} \underset{k_3}{\stackrel{k_1}\leftrightharpoons} 2\mathrm{B} \stackrel{k_2} \longrightarrow \mathrm{ P}$ but you can set $k_3$ to zero if you want. (its called km1 in the code). I've used log time spaced points so that a wide time range can be seen. The x axis is log in the plots.

I have used a built in solver, odeint(), buy you could write your own based on the Euler or Runge -Kutta methods.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def dX_dt(X, t, const):     # returns rates dEf/dt, only used in odeint()

    k1, km1, k2 = const     # unpack constants needed for calc'n
    A = X[0] 
    B = X[1] 
    dAdt=    -k1*A  +  km1*B*B                 #  rate equations here 
    dBdt=     k1*A  - (km1+k2)*B*B

    f = np.array( [(dAdt), (dBdt)] )           # evaluate equations; 
    return f

C0 = np.zeros(2,dtype=float)
C0[0] = 1.0e-3    # C0= [ A(0), B(0)   ]      # initial concentrations 
numt = 200
t = np.logspace( -9,-4,  numt )    # time as initial power, final power, num points
k1  = 1e7  # 1st order rate const s^(-1)
km1 = 1e6  # 2nd dm^3mol^(-1)s^(-1)
k2  = 1e9  # 2nd  ....
const = [ k1, km1, k2 ]                            # rate constants for calc'n

soln = odeint(dX_dt, C0, t, args = (const,)  )     # solve equations                 
Aval, Bval = soln.T     # extract data

plt.plot(t,Aval,color='red',  label='A')


The figure shows the populations vs (log) time abbc scheme

  • $\begingroup$ I'm not sure what to do with your answer because I took your advice - numerical solution but I don't use python so I cannot comment on your code. I'm going to accept your answer because you suggested numerical solution. $\endgroup$
    – 7E10FC9A
    Jun 28, 2018 at 0:14
  • $\begingroup$ You will find that it is quite easy to convert it into your preferred language. Most now do have the equivalent of the odeint integrator used above. Alternatively download python (along with numpy, matplotlib & scipy ) as all are free to use. $\endgroup$
    – porphyrin
    Jul 1, 2018 at 7:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.