\vdots  & \ddots & \vdots \\[3pt]
        d^{-1}_{1N} \quad &  d_{N2}^{-1}\quad & d_{N3}^{-1} \quad& \quad\dots & 0\\[3pt]
           \end{array}
           \right], \quad {\boldsymbol{F}} = \left[\begin{array}{c}
                                          f_1\\[3pt]f_2\\[3pt]\\[3pt]\vdots\\[3pt] f_N
                                         \end{array}\right]$$ Therefore from  we obtain a non-linear system of equations: $${\boldsymbol{F}} - {\boldsymbol{C}} \langle {\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}}\rangle^{3/2} = 0,
 \label{eq:sys}$$ where ${\boldsymbol{C}} = \mathrm{diag}[c_1 \; c_2 \; \dots \; c_N]$ is the diagonal matrix, and the vector ${\boldsymbol{G}}_0 = [g_{10} \; g_{20} \; \dots \; g_{N0}]^\intercal$ contains the local penetrations in absence of elastic interactions. In used notations we also assume that Macaulay brackets applied to a vector and raising to $n$-th power imply the following operation with vector components: $$\langle {\boldsymbol{V}} \rangle = [\langle {\boldsymbol{v}}_1 \rangle\; \langle {\boldsymbol{v}}_2 \rangle\; \dots \; \langle {\boldsymbol{v}}_N \rangle\;]^\intercal,\quad {\boldsymbol{V}}^n = [v_1^n \; v_2^n \; \dots \; v_N^n]^\intercal.$$ Iterative explicit resolution of this problem (i.e. when to compute the configuration on the current load step, forces from the previous load steps are used to estimate the long-range interaction) in most cases leads to big errors. On the other hand, in implicit integration, the system converges badly if Newton-Raphson method is used, because of the discontinuous nature of the function $\langle x\rangle$ and non-coercive properties. Therefore we regularize it as 