Disclaimer: I haven't read most of the papers that the book refers to, so what you will find below is my own understanding based on very little data (i.e. I may have severely misunderstood some part of the book).
Laplace-Beltrami Operator
In a finite element formulation you typically have a stiffness matrix $W$ and a mass matrix $M$. For a 2-manifold one possible discretisation is the cotangent one:
$$W_{ij} =
\begin{cases}
-\frac{1}{2}(\cot\alpha_{ij}+\cot\beta_{ij}) &\text{for } j\in\mathcal{N}(i),\\
\sum_{k\in\mathcal{N}(i)}C_{ik} &\text{for } j=i, \\
0 &\text{for } j\ne i \land j\notin \mathcal{N}(i),
\end{cases}
$$
where $\alpha_{ij}$ and $\beta_{ij}$ are the angles opposite to the edge $ij$. You can also find the same up to a factor $\frac{1}{2}$ in the appendix of the book you linked (note that in said book they use instead of $W$ and $M$, $M$ and $D$), or a few pages prior to the one you cite. In the appendix of the book that you have the mass matrix $M$ is a lumped (diagonal) matrix, where $M_{ii}$ is the Voronoi area for vertex $i$. Then finally they use $\Delta_S\approx L = -M^{-1}W$ (here $\Delta_S$ describes the Laplace-Beltrami operator on the manifold $S$).
Continuous Harmonic Formulation
As far as I can tell the continuous formulation of the fairing approach with a membrane surface would be something like this:
\begin{align}
\Delta_S x &= 0, \quad x\in S, \\
x &= p, \quad x\in \partial S,
\end{align}
where $S$ is the blue part in the images, $\partial\Omega$ is the boundary (the interface between white and blue), and $p$ are the (fixed) positions at said boundary. You may note that the above is actually a non-linear equation. This is the case because $\Delta_{S}$ depends on the surface itself, and $x$ also depends on the surface itself. Thus I would actually suggest (for the sake of devising a numerical algorithm) to instead consider the time evolution equation:
\begin{alignat}{3}
\partial_t x(t) &= \Delta_{S(t)} x(t), &\quad& x\in S(t), &\quad& t\in(0,\infty) \\
x(t) &= p, &\quad& x\in \partial S(t), &\quad& t\in[0,\infty), \\
S(0) &= S_0.
\end{alignat}
This should be equivalent to the first formulation at the steady state $\partial_t x = 0$.
Discrete Harmonic Formulation
The first continuous formulation can be discretised just as:
$$(C-(I-C)L)X = CP,$$
where $\Delta_{S} \approx L = -M^{-1}W\in\mathbb{R}^{n\times n}$, $X\in\mathbb{R}^{n\times 3}$ are the $n$ vertices' positions that are inside the blue region and also the ones from the boundary (vertices that are from the white region but have a blue neighbour), $P\in\mathbb{R}^{n\times 3}$ are the initial positions. Finally $C = \operatorname{diag}(c)$, where $c \approx 1_{\partial S}$:
$$c_i=\begin{cases} 1, \quad \text{$i$ is a boundary vertex}, \\ 0, \quad \text{$i$ is not a boundary vertex}. \end{cases}$$
$(C-(I-C)L)$ essentially negated $L$ and replaces its rows corresponding to Dirichlet vertices with $e_k$ (where $k$ is the index of the row being replaced). You can fully eliminate the Dirichlet vertices by considering the reduced formulation as in FEM also. Note however that in either case $L$ is actually a function of $X$ (since these cotan angles and areas depend on the position of the vertices), this means that the system is non-linear. A standard approach is to try and use a fixed point iteration, i.e. have an initial guess $X^0$, then compute $L^0$ from it, and solve the system $(C-(I-C)L^0)X^1 = CP$. Repeating this yields the iteration $(C-(I-C)L^k)X^{k+1} = CP$, and you stop when e.g. $\|X^{k+1}-X^k\|<\epsilon$ in some iteration (note that the superscripts here are iteration indices and not powers). I have no idea how stable this is however. You can of course try other approaches for solving nonlinear systems of equations.
For the time-evolution equation we can discretise the time derivative, e.g. using explicit or semi-implicit Euler. In the explicit variant you have:
$$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^k,$$
while in the semi-implicit variant you have:
$$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^{k+1}.$$
You can of course also try a fully-implicit variant, but that will result in a nonlinear system of equations that you must solve. Once you discretise also the spatial domain and incorporate the boundary conditions you get:
\begin{alignat}{2}
x^{k+1} &= (I-C)(I+\tau L^k)x^k + CP &\quad&& \text{(explicit Euler)}\\
(C+(I-C)(I-\tau L^k))x^{k+1} &= CP &\quad&& \text{(semi-implicit Euler)}.
\end{alignat}
The semi-implicit variant should be stable for all $\tau>0$ I believe, but larger step sizes may give poor results.
Polyharmonic Formulation
For the polyharmonic formulation you need additional boundary conditions. What I can suggest are Navier conditions:
\begin{align}
(-\Delta_S)^q x &= 0, \quad x\in S, \\
(-\Delta)^l x &= p^l, \quad x\in \partial S, \, 0\leq l \leq q-1.
\end{align}
Then this gets discretised as
$$\left(C\sum_{l=0}^{q-1}(-L)^{l}+(I-C)(-L)^q\right)X = C\sum_{l=0}^{q-1}P^l.$$
Note that the superscripts here refer to powers, except for the superscript of $P$ which refers to just different matrix indices.
Implementation Details
When you solve the linear systems you do not have to solve simultaneously for each coordinate, you can solve them one at a time. The system matrix remains the same, only the right-hand side changes. You may also want symmetrize the $L$ as described in the book, then you could use a conjugate gradient solver since the matrices should be symmetric positive definite (provided you have at least 1 Dirichlet vertex).
When you have to deal with the $(-L)^j$ powers you can just implement them initially as iterated matrix-vector products. Later on you can compute the explicit form of the matrices in order to reuse cache more.
The choice of the Navier conditions for the polyharmonic case lets you directly use $(-L)^l$ (as seen in the last discretisation) instead of having to derive special modification for the boundary terms.
Also, please post some screenshots if you make this work. It's nice to know if what I mentioned actually works out and see the result of it.
