## 3 Numerical Scheme
The identities (9), (10), and (12) inspire the following quasi-Newton scheme given, for a starting value \(x_{0}\in X_{h}\), through the iterative solution of the following minimization problem:  
\[\begin{array}{l}\text{Find }\big{(}c_{S,LS,h}^{(n)},p_{S,LS,h}^{(n)}\big{)}_{S \in\{A,B,C\}}\in X_{h},\text{ such that the following identity holds true}\\ \\ \big{(}v_{S,h}^{(n)},q_{S,h}^{(n)}\big{)}_{S\in\{A,B,C\}}=\underset{\big{(}v_{S,h},d_{S,h}\big{)}_{S\in\{A,B,C\}}}{\text{argmin}}\quad LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}},\mu;\mathcal{T},n\big{)}\\ \\ \text{where }LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}},\mu;\mathcal{T} \big{)}\text{ is, for all }(v_{S,h},q_{S,h})_{S\in\{A,B,C\}}\in X_{h},\text{ defined as}\\ \\ LS_{A}\big{(}(v_{S,h},q_{S,h})_{S\in\{A,B,C\}},\mu;\mathcal{T} \big{)}:=\\ \\ \big{\|}D_{A}\text{div}q_{A}+\gamma_{A}v_{A,h}+k_{2}v_{C,h}-k_{1}\big{[}\Pi_{0 }(c_{A,h}^{(n-1)})v_{B,h}+\Pi_{0}(c_{B,h}^{(n-1)})\
u_{A,h}-\Pi_{0}(c_{A,h}^{( n-1)}c_{B,h}^{(n-1)})\big{]}+\Pi_{0}f_{A}\big{\|}_{L^{2}(\Omega)}^{2}\\ \\ \quad+\big{\|}D_{B}\text{div}q_{B}+\gamma_{B}v_{B,h}+k_{2}v_{C,h}-k_{1}\big{[} \[PERSON]_{0}(c_{A,h}^{(n-1)})v_{B,h}+\[PERSON]_{0}(c_{B,h}^{(n-1)})v_{A,h}-\[PERSON]_{0}(c_{A,h }^{(n-1)}c_{B,h}^{(n-1)})\big{]}+\[PERSON]_{0}f_{B}\big{\|}_{L^{2}(\Omega)}^{2}\\ \\ \quad+\big{\|}D_{C}\text{div}q_{C}+\gamma_{C}v_{C,h}-k_{2}v_{A,h}+k_{1}\big{[} \Pi_{0}(c_{A,h}^{(n-1)})v_{B,h}+\Pi_{0}(c_{B,h}^{(n-1)})v_{A,h}-\Pi_{0}(c_{A,h }^{(n-1)}c_{B,h}^{(n-1)})\big{]}+\Pi_{0}f_{C}\big{\|}_{L^{2}(\Omega)}^{2}\\ \\ \quad+\sum_{S\in\{A,B,C\}}\bigg{[}\big{\|}v_{S}-g_{S,D}\big{\|}_{L^{2}(\Gamma_ {D})}^{2}+\big{\|}\
abla v_{S}\cdot\
u_{N}-g_{S,N}\big{\|}_{L^{2}(\Gamma_{N})} ^{2}\bigg{]}\\ \\ \text{where the set }A_{S}^{(n-1)}\subset\Omega\text{ is the support of }\big{(}\mu_{S}-c_{S}^{(n-1)}\big{)}_{+},\text{ where the support of a function}\\ \quad f\in L^{2}(\Omega)\text{ is defined as }\text{supp}(f):=\overline{\{x\in\Omega|f(x) \
eq 0\}}.\text{ Hence, }A_{S}^{(n-1)}\text{ is defined as:}\\ \\ A_{S}^{(n-1)}=\text{supp}\bigg{[}\big{(}\mu_{S}-c_{S}^{(n-1)}\big{)}_{+}\bigg{]} \\ \text{furthermore for a subset }B\subset\Omega\text{, the map }\chi_{B}:\Omega\to[0,1]\text{ is given by}\\ \\ \chi_{B}(x):= \begin{cases}1,&\text{for }x\in B,\\ 0,&\text{else.}\end{cases} \]  
As it can directly be seen, the inclusion \(\chi_{A(t-1)}(x)\subseteq\partial_{t\in\mathcal{T}}\big{(}[\mu_{S}-c_{S}]( \bullet)\big{)}_{+}(x)\) holds true. As it is commonly known, cf. [17; 45; 46] Newton and quasi-Newton schemes tend to be unstable. One strategy, as described in [46], is to use stabilized quasi-Newton methods based on a coupling of a quasi-Newton method with a subgradient scheme, on a not linearised problem. Another option, which is used and implemented for the numerical examples, as discussed in the Sections 4 and 5, are homotopy methods, cf. [31; 32], which is a scheme for the solution (13), based on the iterative transformation of a system with known solution into the non-linear system associated to (13). The homotopy methodology is probably the easiest way to stabilize the numerical strategy and increase the robustness. **Remark 2**: _The strategy refered to as quasi-Newtonian scheme is in fact a disturbed variant of a 'normal' quasi-Newtonian scheme. The difference between the 'normal' quasi-Newtonian scheme and the quasi-Newtonian scheme discussed in this paper is the application of the projections \(\Pi_{0}\). The sequence that is generated by the corresponding quasi-Newtonian scheme is denoted by \(y^{(k)}\)._  
Some notations are needed to describe the used transformation. First, let, for \(M\in\mathbb{N}\), the following discretisation of the interval \([0,1]\) be given through \(0=t_{1}\leq\ldots\leq t_{j}\leq\ldots\leq t_{M}=1\).