## 3 Numerical Scheme
Subsequently, a homotopy method is given by the solution, for \(1\leq j\leq N\), by computing the following Algorithm 2:```
Input :Initial discretization \(0=t_{1}\leq\ldots\leq t_{j}\leq\ldots\leq t_{M}\) of the interval \([0,1]\), system parameters \(D_{A}\), \(D_{B}\), \(D_{\mathbb{C}}\), \(k_{1}\), \(k_{2}\) and set \(j=0\), set initial value \(x^{(0)}\in X^{h}\), \(0<\varepsilon\). Output :Approximate solution to (8) and set \((c^{(n)}_{LS,S,h},p^{(n)}_{LS,S,h})_{S\in\{A,B,C\}}=(c^{(n,1)}_{LS,S,h},p^{(n,1 )}_{LS,S,h})_{S\in\{A,B,C\}}\). The homotopy algorithm is given through the following steps. S1: For \(n=1,2,\ldots\) solve the following problem. Find \(\big{(}(c^{(n,j)}_{LS,S,h},p^{(n,j)}_{LS,S,h})\big{)}_{S\in\{A,B,C\}}\in X_{h}\) s.t. \[(c^{(n,j)}_{LS,S,h},p^{(n,j)}_{LS,S,h})_{S\in\{A,B,C\}}=\underset{(v_{S,h},q_{ S,h})_{S\in\{A,B,C\}}\in X_{h}}{\operatorname{argmin}} LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}};\mathcal{T},n,j\big{)}\] (14) and break if \[\sum_{S\in\{A,B,C\}}\left[\left\|c^{j,n}_{S,h}-c^{j,n-1}_{S,h}\right\|^{2}_{L ^{2}(\Omega)}+\left\|q^{j,n}_{S,h}-q^{j,n-1}_{S,h}\right\|^{2}_{L^{2}(\Omega)} \right]<\varepsilon\] and for all \((v_{S,h},q_{S,h})_{S\in\{A,B,C\}}\in X_{h}\), the linear functional \(LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}};\mathcal{T},n,j\big{)}\) is given by \[LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}};\mathcal{T},n,j \big{)}\] \[:=(1-t_{j})\sum_{S\in\{A,B,C\}}\left[\left\|v_{S,h}\right\|^{2}_{ L^{2}(\Omega)}+\left\|q_{S,h}\right\|^{2}_{L^{2}(\Omega)}\right]\] \[+t_{j}LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}};\mathcal{T}, n\big{)}.\] S2: If (14) converges for \(n\to\infty\), go to S3 else refine the decomposition \(0\leq t_{0}\leq\ldots\leq t_{M}=1\), set \(j=0\) and go to S1. S3: If \(t_{j}=1\) break the loop, else set \(j=j+1\) and go to S1. ```

**Algorithm 2**Homotopie-Loop

The usage of the homotopy-algorithm, as described above, yields a cost-intensive numerical scheme. Hence, the usage of adaptive schemes for (5) is the best choice at hand to treat the iterate problems. Adaptive schemes in the homotopy steps are possible, but not implemented for this paper. Adaptive FEM is used as the last step of the numerical scheme. As described in [24, 47, 48] the adaptive scheme is given through the following Algorithm 3:

```
Input :Initial regular triangulation \(\mathcal{T}_{0}\), system parameters, set \(j=0\). Output :Fine mesh \(\mathcal{T}_{\ell}\), approximate limit of the solutions to (13), estimation of the error. The adaptive FEM algorithm is given by the following steps: S1: Solve: Take the limit of the solutions iterated problems (13), approximated through Algorithm 2. S2: Estimate: Estimate the error of the current discrete problem on each element \(T\in\mathcal{T}_{j}\) with an estimator \(\eta\). S3: Break: If the triangulation is fine enough or the estimated error is lower than a given Tolerance, then break the condition. S4: Mark: Mark the elements of the triangulation \(T\in\mathcal{T}_{j}\) with the highest estimated local errors. S5: Refine: Refine all marked elements and additonal elements, such that the result \(\mathcal{T}_{j+1}\) is again a regular triangulation. S6: Set \(j=j+1\) and go to S1:At the end of the section, a few remarks have to be made:  
**Remark and Definition 1** (Used setup of the AFEM-Algorithm for this paper). : _Note that:_  
1. _In_ \(S\)2 _of Algorithm_ 3_, an estimator is evaluated on each element_ \(T\in\mathcal{T}_{j}\)_, which has to be defined. Assuming that for the exact solution_ \((u_{LS,S},p_{LS,S})_{S\in\{A,B,C\}}\in X\) _of (_3_), the equality_ \[LS\big{(}(u_{LS,S},p_{LS,S})_{S\in\{A,B,C\}}\big{)}=0\] (15) _holds true.