## 3 Numerical Algorithm of Solving the Problem  
To solve the problem (7)-(15), we use finite difference method [47,48]. In the area \(D=\left\{0\leq x<\infty,\ 0\leq t\leq T\right\}\), we consider the net  
\[\omega_{h\tau}=\left\{\left(x_{k},\,t_{j}\right),\ x_{k}=kh,\ k=0,1,...,\ \ t_{j}=j\tau,\ j=0,1,...,J,\ \ \ \tau=T/J\right\},\]  
where, \(J\) is the number of nodes of time. Instead of functions \(c^{(i)}\left(t,\,x\right)\), \(\rho_{a}^{(i)}\left(t,\,x\right)\), \(\rho_{p}^{(i)}\left(t,\,x\right)\), we consider grid functions, whose values at the nodes \(\left(x_{k},\ t_{j}\right)\) are denoted by \(\left(c^{(i)}\right)_{k^{\prime}}^{j}\left(\rho_{a}^{(i)}\right)_{k^{\prime}} ^{j}\left(\rho_{p}^{(i)}\right)_{k^{\prime}}^{j}\), respectively. The balance Equation (7) are approximated by the equations  
\[ \begin{array}{l}m_{k}^{j}\frac{\left(c^{(i)}\right)_{k}^{j+1}-\left(c^{(i)} \right)_{k}^{j}}{\tau}+v\frac{\left(c^{(i)}\right)_{k}^{j+1}-\left(c^{(i)} \right)_{k-1}^{j+1}}{\hbar}+\frac{\left(\rho_{a}^{(i)}\right)_{k}^{j+1}-\left( \rho_{a}^{(i)}\right)_{k}^{j}}{\tau}+\\ +\frac{\left(\rho_{p}^{(i)}\right)_{k}^{j+1}-\left(\rho_{p}^{(i)}\right)_{k}^{ j}}{\tau}=D^{(i)}\frac{\left(c^{(i)}\right)_{k-1}^{j+1}-2\left(c^{(i)} \right)_{k}^{j+1}+\left(c^{(i)}\right)_{k+1}^{j+1}}{\hbar^{2}}.\end{array}  \tag{16}\]  
Combining the Equations (8) with (9), (10), \(\left|\
abla p\right|\) can be written in grid form as follows  
\[\left|\
abla p\right|_{k}^{j+1}=\frac{v_{0}\left(1-m_{0}+\left(\left(\rho_{a}^ {(1)}\right)_{k}^{j}+\left(\rho_{a}^{(2)}\right)_{k}^{j}+\left(\rho_{n}^{(1)} \right)_{k}^{j}+\left(\rho_{n}^{(2)}\right)_{k}^{j}\right)\right)^{2}}{k_{0} \left(m_{0}-\left(\left(\rho_{a}^{(1)}\right)_{k}^{j}+\left(\rho_{a}^{(2)} \right)_{k}^{j}+\left(\rho_{n}^{(1)}\right)_{k}^{j}+\left(\rho_{n}^{(2)} \right)_{k}^{j}\right)\right)^{3}}. \tag{17}\]  
The difference scheme for the Equation (11) is  
\[\frac{\left(\rho_{a}^{(i)}\right)_{k}^{j+1}-\left(\rho_{a}^{(i)}\right)_{k}^{j} }{\tau}=\frac{\beta_{a}^{(i)}}{1+\gamma^{(i)}\left|\
abla p\right|_{k}^{j+1}} \left(c^{(i)}\right)_{k}^{j}-\beta_{a}^{(i)}\frac{\left(\rho_{a}^{(i)}\right)_{ k}^{j}\left(1+\omega^{(i)}\left|\
abla p\right|_{k}^{j+1}\right)}{\rho_{a0}^{(i)}}c_{0 }^{(i)}. \tag{18}\]For Equations (12)-(14), we have  
\[\frac{\left(\rho_{p}^{(i)}\right)_{k}^{j+1}-\left(\rho_{p}^{(i)}\right)_{k}^{j}}{ \tau}=\rho_{p}^{(i)}\phi_{i}\left(\left(\rho_{p}^{(1)}\right)_{k}^{j},\ \left(\rho_{p}^{(2)}\right)_{k}^{j}\right)\left(c^{(i)}\right)_{k}^{j}, \tag{19}\]  
\[\phi_{1}\left(\left(\rho_{p}^{(1)}\right)_{k}^{j},\ \left(\rho_{p}^{(2)} \right)_{k}^{j}\right)=\left\{ \begin{array}{ll}1&\text{for}\ \ 0<\left(\rho_{p}^{(1)} \right)_{k}^{j}\leq\rho_{p1}^{(1)},\\ \rho_{p1}^{(1)}\Big{/}\left(\rho_{p}^{(1)}\right)_{k}^{j}&\text{for}\ \ \rho_{p1}^{(1)}< \left(\rho_{p}^{(1)}\right)_{k}^{j},\ \left(\rho_{p}^{(1)}\right)_{k}^{j}+ \left(\rho_{n}^{(2)}\right)_{k}^{j}<\rho_{p0},\\ 0&\text{for}\ \ \left(\rho_{p}^{(1)}\right)_{k}^{j}+\left(\rho_{p}^{(2)} \right)_{k}^{j}=\rho_{p0},\end{array} \right. \tag{20}\]  
\[\phi_{2}\left(\left(\rho_{n}^{(1)}\right)_{k}^{j},\ \left(\rho_{n}^{(2)} \right)_{k}^{j}\right)=\left\{ \begin{array}{ll}1&\text{for}\ \ 0<\left(\rho_{p}^{(2)} \right)_{k}^{j}\leq\rho_{p1}^{(1)},\\ \rho_{p1}^{(2)}\Big{/}\left(\rho_{p}^{(2)}\right)_{k}^{j},\ \text{for},\rho_{p1}^{(2)}< \left(\rho_{p}^{(2)}\right)_{k}^{j}<\rho_{p0}^{(2)},\left(\rho_{p}^{(1)} \right)_{k}^{j}+\left(\rho_{p}^{(2)}\right)_{k}^{j}<\rho_{p0},\\ 0&\text{for}\ \ \left(\rho_{p}^{(2)}\right)_{k}^{j}=\rho_{p0}^{(2)}\ \text{or}\ \ \left(\rho_{p}^{(1)} \right)_{k}^{j}+\left(\rho_{p}^{(2)}\right)_{k}^{j}=\rho_{p0}.\end{array} \right.