A Proof of Accuracy for the MPRK Scheme (3.25)
In this appendix, we show that the MPRK scheme (3.25) is third-order accurate for the system of ODEs (3.24).
According to (3.25a) and (3.25b), \(c_{k,i}^{(1)}\) can be expanded as
$$\begin{aligned} c_{k,i}^{(1)} = c_{k,i}^n + \Delta t \beta _{10} S_{k,i} + \Delta t^2 \beta _{10}^2 Q_{k,i} + \Delta t^3 \beta _{10}^3 R_{k,i} + {\mathcal {O}}(\Delta t^4), \end{aligned}$$
(A.1)
with
$$\begin{aligned} \begin{aligned} S_{k,i}&:= F_{k,i}^n + P_{k,i}^n - D_{k,i}^n , \\ Q_{k,i}&:= \sum _j p_{k,ij}^n \frac{F_{k,j}^n + P_{k,j}^n - D_{k,j}^n}{c_{k,j}^n} - \sum _j d_{k,ij}^n \frac{F_{k,i}^n + P_{k,i}^n - D_{k,i}^n}{c_{k,i}^n} , \\ R_{k,i}&:= \left( \sum _{j,r} \frac{p_{k,ij}^n ( p_{k,jr}^n \frac{F_{k,r}^n + P_{k,r}^n - D_{k,r}^n}{c_{k,r}^n} - d_{k,jr}^n \frac{F_{k,j}^n + P_{k,j}^n - D_{k,j}^n}{c_{k,j}^n} )}{c_{k,j}^n} \right. \\&\quad \left. - \sum _{j,r} \frac{d_{k,ij}^n ( p_{k,ir}^n \frac{F_{k,r}^n + P_{k,r}^n - D_{k,r}^n}{c_{k,r}^n}- d_{k,ir}^n \frac{F_{k,i}^n + P_{k,i}^n - D_{k,i}^n}{c_{k,i}^n} )}{c_{k,i}^n} \right) . \end{aligned} \end{aligned}$$
Thus
$$\begin{aligned} \begin{aligned}&\frac{\rho _{k,i}}{c_{k,i}^n} = n_1 \frac{c_{k,i}^{(1)}}{c_{k,i}^n} + n_2 ( \frac{c_{k,i}^{(1)}}{c_{k,i}^n} )^2 \\&\quad = 1 + (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21} - x) \Delta t \frac{S_{k,i}}{c_{k,i}^n} + (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21} - x) \beta _{10} Q_{k,i} \Delta t^2 \\&\qquad + (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21} - x - \beta _{10})\beta _{10} \left( \frac{S_{k,i}}{c_{k,i}^n} \right) ^2 \Delta t^2 + {\mathcal {O}}(\Delta t^3). \end{aligned} \end{aligned}$$
(A.2)
Recall from (2.24) that \(\phi ^{(l)}\) with \(l=1,2\) and \(\phi =p_{k,ij}, d_{k,ij}, F_{k,i}\) can be expanded as
$$\begin{aligned} \phi ^{(l)} = \phi (\mathbf {c}^{(l)}) = \phi ^n + \frac{\partial \phi ^n }{\partial \mathbf {c}} (\mathbf {c}^{(l)} - \mathbf {c}^n) + \frac{1}{2}(\mathbf {c}^{(l)} - \mathbf {c}^n)^T \mathbf {H}_{\phi }^n (\mathbf {c}^{(l)} - \mathbf {c}^n) + {\mathcal {O}}(\Delta t^3). \end{aligned}$$
(A.3)
Taking \(l=1\) and substituting (A.1) into the above equation, we have
$$\begin{aligned} \phi ^{(1)} = \phi ^n + \Delta t \beta _{10} \frac{\partial \phi ^n }{\partial \mathbf {c}}\mathbf {S} + \Delta t^2 \beta _{10}^2 \frac{\partial \phi ^n }{\partial \mathbf {c}} \mathbf {Q} + \frac{1}{2}\beta _{10}^2 \Delta t^2 \mathbf {S}^T \mathbf {H}_{\phi }^n \mathbf {S} + {\mathcal {O}}(\Delta t^3) \end{aligned}$$
(A.4)
with \(\phi =p_{k,ij}, d_{k,ij}, F_{k,i}\). It follows from this, (3.25d) and (A.2) that
$$\begin{aligned} \begin{aligned} c_{k,i}^{(2)}&= \alpha _{20}c_{k,i}^{(0)} + \alpha _{21}c_{k,i}^{(1)} + \Delta t \beta _{20} F_{k,i}(\mathbf {c}^{(0)})+ \Delta t \beta _{21} F_{k,i}(\mathbf {c}^{(1)}) \\&\quad +\Delta t \sum _j \left( (\beta _{20} p_{k,ij}^{(0)} + \beta _{21} p_{k,ij}^{(1)} ) \frac{c_{k,j}^{(2)}}{\rho _{k,j}} - ( \beta _{20}d_{k,ij}^{(0)} + \beta _{21}d_{k,ij}^{(1)} ) \frac{c_{k,i}^{(2)}}{\rho _{k,i}} \right) \\&= c_{k,i}^n + (\alpha _{21} \beta _{10} + \beta _{20} + \beta _{21} ) S_{k,i} \Delta t + {\mathcal {O}}(\Delta t^2) \end{aligned} \end{aligned}$$
and
$$\begin{aligned} \frac{c_{k,i}^{(2)}}{\rho _{k,i}} = c_{k,i}^n + x S_{k,i} \Delta t + {\mathcal {O}}(\Delta t^2) . \end{aligned}$$
The above two equations further yield
$$\begin{aligned} c_{k,i}^{(2)}= & {} c_{k,i}^n + (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) \Delta t S_{k,i} + \beta _{21} \beta _{10} \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}}\mathbf {S} \\&+ \,\Delta t^2[ \alpha _{21}\beta _{10}^2 + (\beta _{20}+\beta _{21}) x]Q_{k,i} + {\mathcal {O}}(\Delta t^3) \end{aligned}$$
and
$$\begin{aligned} \frac{c_{k,i}^{(2)}}{\rho _{k,i}} = c_{k,i}^n + x S_{k,i} \Delta t + T_{k,i}\Delta t^2 + {\mathcal {O}}(\Delta t^3) \end{aligned}$$
with
$$\begin{aligned} \begin{aligned}&T_{k,i} = \left[ - \beta _{20}\beta _{10} - \beta _{21}\beta _{10} + (\beta _{10} + \beta _{20} + \beta _{21}) x \right] Q_{k,i} + \beta _{21} \beta _{10} \frac{\partial S_{k,i}}{\partial \mathbf {c}} \frac{\mathbf {S}}{c_{k,i}^n} \\&\quad - \left[ (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21} - x - \beta _{10}) \beta _{10} + x( \alpha _{21}\beta _{10} + \beta _{20} + \beta _{21} - x ) \right] \left( \frac{S_{k,i}}{c_{k,i}^n} \right) ^2 . \end{aligned}\qquad \end{aligned}$$
(A.5)
Furthermore, we denote
$$\begin{aligned} X_{k,i} := x \frac{S_{k,i}}{c_{k,i}^n} \end{aligned}$$
and have the expansion of \(c_{k,i}^{(2)}\):
$$\begin{aligned}&c_{k,i}^{(2)} = \alpha _{20}c_{k,i}^{(0)} + \alpha _{21}c_{k,i}^{(1)} + \Delta t \beta _{20} F_{k,i}(\mathbf {c}^{(0)})+ \Delta t \beta _{21} F_{k,i}(\mathbf {c}^{(1)}) \nonumber \\&\qquad +\Delta t \sum _j \left( (\beta _{20} p_{k,ij}^{(0)} + \beta _{21} p_{k,ij}^{(1)} ) \frac{c_{k,j}^{(2)}}{\rho _{k,j}} - ( \beta _{20}d_{k,ij}^{(0)} + \beta _{21}d_{k,ij}^{(1)} ) \frac{c_{k,i}^{(2)}}{\rho _{k,i}} \right) \nonumber \\&\quad = c_{k,i}^n + \Delta t \alpha _{21} \beta _{10} S_{k,i} + \Delta t^2 \alpha _{21} \beta _{10}^2 Q_{k,i} + \Delta t^3 \alpha _{21} \beta _{10}^3 R_{k,i}\nonumber \\&\qquad + \Delta t( \beta _{20} + \beta _{21})F_{k,i}^n +\Delta t^2 \beta _{10}\beta _{21} \frac{\partial F_{k,i}^n }{\partial \mathbf {c}}\mathbf {S}\nonumber \\&\qquad + \Delta t^3 \beta _{10}^2 \beta _{21} \frac{\partial F_{k,i}^n }{\partial \mathbf {c}} \mathbf {Q} + \frac{1}{2} \beta _{10}^2 \beta _{21} \Delta t^3 \mathbf {S}^T \mathbf {H}_{F_{k,i}}^n \mathbf {S} \nonumber \\&\qquad + \beta _{20} \Delta t \sum _j p_{k,ij}^n(1 + X_{k,j} \Delta t + T_{k,j} \Delta t^2 ) - \beta _{20} \Delta t \sum _j d_{k,ij}^n(1 + X_{k,i} \Delta t + T_{k,i} \Delta t^2 ) \nonumber \\&\qquad + \beta _{21} \Delta t \sum _j p_{k,ij}^n + \beta _{21} \beta _{10} \Delta t^2 \sum _j \frac{\partial p_{k,ij}^n }{\partial \mathbf {c}} \mathbf {S} + \beta _{21} \beta _{10}^2 \Delta t^3 \sum _{j}\frac{\partial p_{k,ij}^n }{\partial \mathbf {c}}\mathbf {Q}\nonumber \\&\qquad + \frac{1}{2} \beta _{21} \beta _{10}^2 \Delta t^3 \sum _j\mathbf {S}^T \mathbf {H}_{p_{k,ij}}^n \mathbf {S} \nonumber \\&\qquad + \beta _{21}\Delta t^2 \sum _j p_{k,ij}^n X_{k,j} + \beta _{21}\beta _{10}\Delta t^3 \sum _j \frac{\partial p_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} X_{k,j} + \beta _{21}\Delta t^3 \sum _j p_{k,ij}^n T_{k,j} \nonumber \\&\qquad - \beta _{21} \Delta t \sum _j d_{k,ij}^n + \beta _{21} \beta _{10} \Delta t^2 \sum _j \frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} - \beta _{21} \beta _{10}^2 \Delta t^3 \sum _{j}\frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {Q}\nonumber \\&\qquad - \frac{1}{2} \beta _{21} \beta _{10}^2 \Delta t^3 \sum _j\mathbf {S}^T \mathbf {H}_{d_{k,ij}}^n \mathbf {S} \nonumber \\&\qquad - \beta _{21}\Delta t^2 \sum _j d_{k,ij}^n X_{k,i} + \beta _{21}\beta _{10}\Delta t^3 \sum _j \frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} X_{k,i} + \beta _{21}\Delta t^3 \sum _j d_{k,ij}^n T_{k,i} + {\mathcal {O}}(\Delta t^4) \nonumber \\&\quad = c_{k,i}^n + (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) \Delta t S_{k,i} + \beta _{21} \beta _{10} \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}}\mathbf {S}\nonumber \\&\qquad + [ \alpha _{21}\beta _{10}^2 + (\beta _{20}+\beta _{21}) x] \Delta t^2 Q_{k,i} \nonumber \\&\qquad + \Delta t^3 \alpha _{21} \beta _{10}^3 R_{k,i} + \Delta t^3 (\beta _{20}+\beta _{21}) ( \sum _j p_{k,ij}^n T_{k,j} - \sum _j d_{k,ij}^n T_{k,i} ) + \beta _{21} \beta _{10}^2 \Delta t^3 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {Q} \nonumber \\&\qquad + \frac{1}{2} \beta _{21} \beta _{10}^2 \Delta t^3 \mathbf {S}^T (\mathbf {H}_{F_{k,i}}+\sum _j\mathbf {H}_{p_{k,ij}}^n-\sum _j\mathbf {H}_{d_{k,ij}}^n) \mathbf {S} \nonumber \\&\qquad + \beta _{21}\beta _{10}x \Delta t^3 \sum _j \left( \frac{\partial p_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,j}}{c_{k,j}^n} - \frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,i}}{c_{k,i}^n} \right) + {\mathcal {O}}(\Delta t^4). \end{aligned}$$
(A.6)
Next we expand \(a_{k,i}\) according to (3.25e) and (3.25f). Note that
$$\begin{aligned} \mu _{k,i} = c_{k,i}^n (\frac{c_{k,i}^{(1)}}{c_{k,i}^n})^s = c_{k,i}^n + s \beta _{10} S_{k,i} \Delta t + {\mathcal {O}}(\Delta t^2) \end{aligned}$$
and thereby
$$\begin{aligned} \begin{aligned} a_{k,i}&= \eta _{1} c_{k,i}^n + \eta _{2} c_{k,i}^{(1)} + \Delta t \eta _{3}(\eta _{1} + \eta _{2} ) F_{k,i}(\mathbf {c}^{(0)}) + \Delta t \eta _{4}(\eta _{1} + \eta _{2} ) F_{k,i}(\mathbf {c}^{(1)})\\&\quad + \Delta t \sum _j \bigg ( \eta _{3}p_{k,ij}^n + \eta _{4}p_{k,ij}^{(1)}\bigg ) \frac{a_{k,j}}{\mu _{k,j}} -\Delta t \sum _j \bigg ( \eta _{3}d_{k,ij}^n + \eta _{4}d_{k,ij}^{(1)}\bigg ) \frac{a_{k,i}}{\mu _{k,i}}, \\&= (\eta _{1} + \eta _{2} ) c_{k,i}^n + \Delta t \eta _2 \beta _{10} S_{k,i} + \Delta t(\eta _{1} + \eta _{2} )(\eta _{3} + \eta _{4} ) F_{k,i}^n \\&\quad + (\eta _{1} + \eta _{2} )(\eta _{3} + \eta _{4} )( P_{k,i}^n - D_{k,i}^n) \Delta t + {\mathcal {O}}(\Delta t^2) \\&= (\eta _{1} + \eta _{2} ) c_{k,i}^n + [ \eta _2 \beta _{10} + (\eta _{1} + \eta _{2} )(\eta _{3} + \eta _{4} )] S_{k,i} \Delta t + {\mathcal {O}}(\Delta t^2) , \\ \frac{a_{k,i}}{\mu _{k,i}} =&\eta _{1} + \eta _{2} + [ \eta _2 \beta _{10} + (\eta _1+\eta _2)(\eta _3+\eta _4 -s \beta _{10} ) ] \frac{S_{k,i}}{c_{i,k}^n} \Delta t + {\mathcal {O}}(\Delta t^2). \end{aligned} \end{aligned}$$
With these, we further expand \(\frac{a_{k,i}}{c_{k,i}^n}\) as
$$\begin{aligned} \begin{aligned} \frac{a_{k,i}}{c_{k,i}^n}&= \eta _{1} + \eta _{2} + [ \eta _2 \beta _{10} + (\eta _1+\eta _2)(\eta _3+\eta _4) ] \frac{S_{k,i}}{c_{i,k}^n} \Delta t \\&\qquad + \{ \eta _2 \beta _{10}^2 + [ \eta _2 \beta _{10} + (\eta _1+\eta _2)(\eta _3+\eta _4 -s \beta _{10} ) ](\eta _{3} + \eta _{4}) \} Q_{k,i} \Delta t^2 \\&\qquad + \eta _4 \beta _{10}( \eta _{1} + \eta _{2} )\frac{ \partial S_{k,i}}{\partial \mathbf {c}} \frac{\mathbf {S}}{c_{k,i}^n} \Delta t^2 + {\mathcal {O}}(\Delta t^3) \\&\quad = r_1 + r_2 \frac{S_{k,i}}{c_{i,k}^n} \Delta t + r_3 Q_{k,i} \Delta t^2 + r_4 \frac{ \partial S_{k,i}}{\partial \mathbf {c}} \frac{\mathbf {S}}{c_{k,i}^n} \Delta t^2 + {\mathcal {O}}(\Delta t^3) , \end{aligned} \end{aligned}$$
(A.7)
where \(r_i,i=1,2,3,4\) are defined below (3.17). It follows from this and (3.13) that
$$\begin{aligned} \begin{aligned} \sigma _{k,i}&= a_{k,i} + z c_{k,i}^{(0)} \frac{c_{k,i}^{(2)}}{\rho _{k,i}} \\&= 1 + (1-y) S_{k,i} \Delta t + \frac{ (\alpha _{31} + \alpha _{32}\alpha _{21})\beta _{10}^3 }{ \beta _{30} + \beta _{31} + \beta _{32} } Q_{k,i} \Delta t^2 + \frac{ \alpha _{32}( \beta _{20} + \beta _{21} ) }{ \beta _{30} + \beta _{31} + \beta _{32} } T_{k,i} \Delta t^2 \\&\quad - y(1-y) c_{k,i}^n \left( \frac{S_{k,i}}{c_{k,i}^n} \right) ^2 \Delta t^2 + \frac{1}{2} \frac{\partial S_{k,i}}{\partial \mathbf {c}} \mathbf {S} \Delta t^2 + {\mathcal {O}}(\Delta t^3). \end{aligned} \end{aligned}$$
(A.8)
Moreover, we substitute expansion (A.6) into (A.3) with \(l=2\) and obtain
$$\begin{aligned} \begin{aligned} \phi ^{(2)}&= \phi ^n + \Delta t (\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) \frac{\partial \phi ^n }{\partial \mathbf {c}}\mathbf {S} + \beta _{21} \beta _{10} \Delta t^2 \sum _{m,r}\frac{\partial \phi ^n }{\partial c_{m,r}} \frac{\partial S_{m,r} }{\partial \mathbf {c}}\mathbf {S}\\&\quad + \Delta t^2 [ \alpha _{21}\beta _{10}^2 + (\beta _{20}+\beta _{21}) x] \frac{\partial \phi ^n }{\partial \mathbf {c}} \mathbf {Q} \\&\quad + \frac{1}{2}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21})^2 \Delta t^2 \mathbf {S}^T \mathbf {H}_{\phi }^n \mathbf {S} + {\mathcal {O}}(\Delta t^3) \end{aligned} \end{aligned}$$
(A.9)
with \(\phi =p_{k,ij}, d_{k,ij}, F_{k,i}\). From this and (A.4) it follows that
$$\begin{aligned} \begin{aligned}&\beta _{30} \phi ^n + \beta _{31} \phi ^{(1)} + \beta _{32} \phi ^{(2)} \\&\quad = (\beta _{30} + \beta _{31} + \beta _{32}) \phi ^n + \Delta t [\beta _{31} \beta _{10} + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) ] \frac{\partial \phi ^n }{\partial \mathbf {c}} \mathbf {S} \\&\qquad + \beta _{32} \beta _{21} \beta _{10} \Delta t^2 \sum _{m,r}\frac{\partial \phi ^n }{\partial c_{m,r}} \frac{\partial S_{m,r} }{\partial \mathbf {c}}\mathbf {S} \\&\qquad + \Delta t^2 [ \beta _{31} \beta _{10}^2 + \beta _{32}\alpha _{21}\beta _{10}^2 + \beta _{32} (\beta _{20}+\beta _{21}) x] \frac{\partial \phi ^n }{\partial \mathbf {c}} \mathbf {Q}\\&\qquad + \frac{1}{2}[ \beta _{31}\beta _{10}^2 + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21})^2] \Delta t^2 \mathbf {S}^T \mathbf {H}_{\phi }^n \mathbf {S} + {\mathcal {O}}(\Delta t^3). \end{aligned} \end{aligned}$$
(A.10)
In addition, with expansions (A.1) and (A.6), we have
$$\begin{aligned}&\alpha _{30} c_{k,i}^n + \alpha _{31} c_{k,i}^{(1)} + \alpha _{32} c_{k,i}^{(2)}\nonumber \\&\quad = c_i^n + \Delta t [\alpha _{31} \beta _{10} + \alpha _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21})] S_{k,i} + \alpha _{32} \beta _{21} \beta _{10} \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}}\mathbf {S} \nonumber \\&\qquad + [ \alpha _{31}\beta _{10}^2 + \alpha _{32} \alpha _{21}\beta _{10}^2 + \alpha _{32}(\beta _{20}+\beta _{21}) x] \Delta t^2 Q_{k,i} + ( \alpha _{31} + \alpha _{32} \alpha _{21})\beta _{10}^3 \Delta t^3 R_{k,i} \nonumber \\&\qquad + \Delta t^3 \alpha _{32}(\beta _{20}+\beta _{21}) ( \sum _j p_{k,ij}^n T_{k,j} - \sum _j d_{k,ij}^n T_{k,i} ) + \alpha _{32} \beta _{21} \beta _{10}^2 \Delta t^3 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {Q} \nonumber \\&\qquad + \frac{1}{2} \alpha _{32}\beta _{21} \beta _{10}^2 \Delta t^3 \mathbf {S}^T (\mathbf {H}_{F_{k,i}}+\sum _j \mathbf {H}_{p_{k,ij}}^n - \sum _j \mathbf {H}_{d_{k,ij}}^n) \mathbf {S} \nonumber \\&\qquad + \alpha _{32}\beta _{21}\beta _{10}x \Delta t^3 \sum _j \left( \frac{\partial p_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,j}}{c_{k,j}^n} - \frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,i}}{c_{k,i}^n} \right) + {\mathcal {O}}(\Delta t^4). \end{aligned}$$
(A.11)
Using this and (A.10), we deduce from (3.25f) that
$$\begin{aligned} \begin{aligned} c_{k,i}^{n+1}&= \alpha _{30}c_{k,i}^{(0)} + \alpha _{31}c_{k,i}^{(1)} + \alpha _{32}c_{k,i}^{(2)} + \Delta t \beta _{30} F_{k,i}(\mathbf {c}^{(0)})+ \Delta t \beta _{31} F_{k,i}(\mathbf {c}^{(1)}) + \Delta t \beta _{32} F_{k,i}(\mathbf {c}^{(2)})\\&\qquad +\Delta t \sum _j \left( \bigg (\beta _{30} p_{k,ij}^{(0)} + \beta _{31} p_{k,ij}^{(1)} + \beta _{32} p_{k,ij}^{(2)} \bigg ) \frac{c_{k,j}^{n+1}}{\sigma _{k,j}}\right. \\&\left. \qquad - \bigg ( \beta _{30}d_{k,ij}^{(0)} + \beta _{31}d_{k,ij}^{(1)} + \beta _{32}d_{k,ij}^{(2)} \bigg ) \frac{c_{k,i}^{n+1}}{\sigma _{k,i}} \right) \\&\quad = c_{k,i}^n + \Delta t [\alpha _{31} \beta _{10} + \alpha _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) + (\beta _{30} + \beta _{31} + \beta _{32})]S_{k,i} + {\mathcal {O}}(\Delta t^2) \\&\quad = c_{k,i}^n + \Delta t S_{k,i} + {\mathcal {O}}(\Delta t^2) \end{aligned} \end{aligned}$$
and
$$\begin{aligned} \frac{c_{k,i}^{n+1}}{\sigma _{k,i}} = 1 + y \frac{ S_{k,i} }{c_{k,i}^n} \Delta t + {\mathcal {O}}(\Delta t^2) . \end{aligned}$$
The above two equations yield
$$\begin{aligned} \begin{aligned} c_{k,i}^{n+1}&= c_i^n + \Delta t S_{k,i} \\&\qquad +[\alpha _{32} \beta _{21} \beta _{10} + \beta _{31} \beta _{10} + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) ] \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {S} + {\mathcal {O}}(\Delta t^3) \\&= c_i^n + \Delta t S_{k,i} + \frac{1}{2}\frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {S} \Delta t^2 + {\mathcal {O}}(\Delta t^3) \end{aligned} \end{aligned}$$
(A.12)
and
$$\begin{aligned} \begin{aligned} \frac{c_{k,i}^{n+1}}{\sigma _{k,i}} = 1 + y \frac{ S_{k,i} }{c_{k,i}^n} \Delta t - \frac{ (\alpha _{31} + \alpha _{32}\alpha _{21})\beta _{10}^3 }{ \beta _{30} + \beta _{31} + \beta _{32} } Q_{k,i} \Delta t^2 - \frac{ \alpha _{32}( \beta _{20} + \beta _{21} ) }{ \beta _{30} + \beta _{31} + \beta _{32} } T_{k,i} \Delta t^2 + {\mathcal {O}}(\Delta t^3) . \end{aligned} \end{aligned}$$
Denote
$$\begin{aligned} N_{k,i} = - \frac{ (\alpha _{31} + \alpha _{32}\alpha _{21})\beta _{10}^3 }{ \beta _{30} + \beta _{31} + \beta _{32} } Q_{k,i} - \frac{ \alpha _{32}( \beta _{20} + \beta _{21} ) }{ \beta _{30} + \beta _{31} + \beta _{32} } T_{k,i} \end{aligned}$$
and thus
$$\begin{aligned} \frac{c_{k,i}^{n+1}}{\sigma _{k,i}} = 1 + y \frac{ S_{k,i} }{c_{k,i}^n} \Delta t + N_{k,i}\Delta t^2 + {\mathcal {O}}(\Delta t^3) . \end{aligned}$$
(A.13)
Then we have the expansion of \(c_{k,i}^{n+1}\):
$$\begin{aligned}&c_{k,i}^{n+1} = \alpha _{30}c_{k,i}^{(0)} + \alpha _{31}c_{k,i}^{(1)} + \alpha _{32}c_{k,i}^{(2)} + \Delta t \beta _{30} F_{k,i}(\mathbf {c}^{(0)})+ \Delta t \beta _{31} F_{k,i}(\mathbf {c}^{(1)}) + \Delta t \beta _{32} F_{k,i}(\mathbf {c}^{(2)})\\&\qquad +\Delta t \sum _j \left( \bigg (\beta _{30} p_{k,ij}^{(0)} + \beta _{31} p_{k,ij}^{(1)} + \beta _{32} p_{k,ij}^{(2)} \bigg ) \frac{c_{k,j}^{n+1}}{\sigma _{k,j}}\right. \\&\left. \qquad - \bigg ( \beta _{30}d_{k,ij}^{(0)} + \beta _{31}d_{k,ij}^{(1)} + \beta _{32}d_{k,ij}^{(2)} \bigg ) \frac{c_{k,i}^{n+1}}{\sigma _{k,i}} \right) \\&\quad = c_i^n + \Delta t [\alpha _{31} \beta _{10} + \alpha _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) + (\beta _{30} + \beta _{31} + \beta _{32})] S_{k,i}\\&\qquad + [\alpha _{32} \beta _{21} \beta _{10} + \beta _{31} \beta _{10} + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) ] \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {S} \\&\qquad + \beta _{32} \beta _{21} \beta _{10} \Delta t^3 \sum _{m,r}\frac{\partial S_{k,i} }{\partial c_{m,r}} \frac{\partial S_{m,r} }{\partial \mathbf {c}}\mathbf {S}\\&\qquad + [ \alpha _{31}\beta _{10}^2 + \alpha _{32} \alpha _{21}\beta _{10}^2 + \alpha _{32}(\beta _{20}+\beta _{21}) x + (\beta _{30} + \beta _{31} + \beta _{32})y] \Delta t^2 Q_{k,i}\\&\qquad +\Delta t^3 \sum _j p_{k,ij}^n [ ( \alpha _{31} + \alpha _{32} \alpha _{21})\beta _{10}^3 \sum _r \frac{\bigg ( p_{k,jr}^n \frac{F_{k,r}^n + P_{k,r}^n - D_{k,r}^n}{c_{k,r}^n} - d_{k,jr}^n \frac{F_{k,j}^n + P_{k,j}^n - D_{k,j}^n}{c_{k,j}^n} \bigg )}{c_{k,j}^n} \\&\qquad + \alpha _{32}(\beta _{20}+\beta _{21})T_{k,j} + (\beta _{30} + \beta _{31} + \beta _{32}) N_{k,j} ] \\&\qquad -\Delta t^3 \sum _j d_{ij}^n [ ( \alpha _{31} + \alpha _{32} \alpha _{21})\beta _{10}^3 \sum _r \frac{\bigg ( p_{k,ir}^n \frac{F_{k,r}^n + P_{k,r}^n - D_{k,r}^n}{c_{k,r}^n} - d_{k,ir}^n \frac{F_{k,i}^n + P_{k,i}^n - D_{k,i}^n}{c_{k,i}^n} \bigg )}{c_{k,i}^n} \\&\qquad + \alpha _{32}(\beta _{20}+\beta _{21})T_{k,i} + (\beta _{30} + \beta _{31} + \beta _{32}) N_{k,i} ] \\&\qquad + \frac{1}{2} [ \beta _{31}\beta _{10}^2 + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21})^2 \\&\qquad + \alpha _{32}\beta _{21} \beta _{10}^2] \Delta t^3 \mathbf {S}^T (\mathbf {H}_{F_{k,i}}^n + \sum _j \mathbf {H}_{p_{k,ij}}^n - \sum _j\mathbf {H}_{d_{k,ij}}^n) \mathbf {S} \\&\qquad + [ \beta _{31} \beta _{10}^2 + \beta _{32}\alpha _{21}\beta _{10}^2 + \beta _{32} (\beta _{20}+\beta _{21}) x + \alpha _{32}\beta _{21} \beta _{10}^2 ] \Delta t^3 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {Q}\\&\qquad + [ \alpha _{32}\beta _{21}\beta _{10}x + \beta _{31} \beta _{10}y + \beta _{32}(\alpha _{21}\beta _{10} + \beta _{20} + \beta _{21}) y] \Delta t^3 \\&\qquad \sum _j\left( \frac{\partial p_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,j}}{c_{k,j}^n} - \frac{\partial d_{k,ij}^n }{\partial \mathbf {c}}\mathbf {S} \frac{S_{k,i}}{c_{k,i}^n} \right) +{\mathcal {O}}(\Delta t^4 )\\&\quad = c_{k,i}^n + \Delta t S_{k,i} + \frac{1}{2} \Delta t^2 \frac{\partial S_{k,i} }{\partial \mathbf {c}} \mathbf {S} + \frac{1}{6} \Delta t^3 \sum _{m,r}\frac{\partial S_{k,i} }{\partial c_{m,r}} \frac{\partial S_{m,r} }{\partial \mathbf {c}}\mathbf {S}\\&\qquad + \frac{1}{6} \Delta t^3 \mathbf {S}^T (\mathbf {H}_{F_{k,i}}^n + \sum _j \mathbf {H}_{p_{k,ij}}^n - \sum _j\mathbf {H}_{d_{k,ij}}^n) \mathbf {S}+ {\mathcal {O}}( \Delta t^4 ). \end{aligned}$$
In the last equality, the additional condition (2.31c) has been used. The above expansion shows that scheme (3.25) for the system of ODEs (3.24) is third-order accurate.
