The positive-interior MPRK22(1) step

Let $(p_{ij},d_{ij})$ be the rates of a finite conservative PDS Finite conservative production destruction systems with state domain $\mathrm{\Omega }$. Given $y^n\in {\mathbb{R}}^m$ and $h$, if $y^n\mathrm{\notin }{\mathrm{\Omega }}_{+}$, any coordinate of $y^n$ is nonpositive, or $h\le 0$, return an inadmissible-input status before evaluating rates or solving a system. First compute the stage $z$ by invoking Modified Patankar Euler for conservative transfer systems with denominators ${\sigma }_i=y_i^n$:

$z_i=y_i^n+h{\sum }_{j\mathrm{\ne }i}(p_{ij}(y^n)\frac{z_j}{y_j^n}-d_{ij}(y^n)\frac{z_i}{y_i^n})\mathrm{.}$

If the stage solve does not return a unique vector, return stage-failure. If its output $z\mathrm{\notin }{\mathrm{\Omega }}_{+}$, return rate-domain-failure before evaluating rates at $z$. Otherwise form $B(y^n,z)$ with, for $i\mathrm{\ne }j$,

$B_{ij}(y^n,z):=\frac{p_{ij}(y^n)+p_{ij}(z)}{2z_j},B_{jj}(y^n,z):=-{\sum }_{i\mathrm{\ne }j}{B}_{ij}(y^n,z),$

and solve

$(I-hB(y^n,z))y^{n+1}=y^n\mathrm{.}$

Equivalently,

$y_i^{n+1}=y_i^n+\frac{h}{2}{\sum }_{j\mathrm{\ne }i}((p_{ij}(y^n)+p_{ij}(z))\frac{y_j^{n+1}}{z_j}-(d_{ij}(y^n)+d_{ij}(z))\frac{y_i^{n+1}}{z_i})\mathrm{.}$

If the final solve is not uniquely solvable, return update-failure; otherwise return its unique solution $y^{n+1}$. No zero-stage or denominator-regularized second-order variant is defined here.