There's two kinds of sampling here, time sampling and spatial sampling, since the string equation is a PDE. I cannot read German, so I cannot understand your article. I will be presenting the discussion from the book, Physical Audio Signal Processing.
The ideal string PDE is $Ky'' = \epsilon\ddot{y}$,
where $K$ is the string tension, $\epsilon$ is the mass density and $y$ is the string displacement.
Now if you use the FDTD method to solve this differential equation, the first and second order partial derivatives are
$$y'(t,x) = \frac{y(t,x+X) - y(t,x)}{X} \\
y''(t,x) = \frac {y(t,x+X) - 2y(t,x) + y(t, x-X)}{X^2} \\
\dot{y}(t,x) = \frac{y(t+T,x) - y(t,x)}{T} \\
\ddot{y}(t,x) = \frac{y(t+T,x) - 2y(t,x) + y(t-T,x)}{T^2}
$$
Here, $T$ is the time-domain sampling rate, and it can be chosen according to Keegs' answer, i.e., twice the frequency of the highest frequency string harmonic. To model all audible string modes, choose $T = \frac{1}{44,100}$s. $X = cT$ is the spatial sampling rate. Here $c = \sqrt{K/\epsilon}$ is the speed of sound.
Now, using $t = nT$ and $x = mX$, and using $c^2 y'' = \ddot{y}$, we have
$$c^2 \frac{y(nT, mX+X) - 2y(nT, mX) + y(nT, mX-X)}{X^2} = \frac{y(nT+T,mX) - 2y(nT,mX)+y(nT-T,mX)}{T^2} \\
c^2\frac{y(n,m+1) - 2y(n,m) + y(n,m-1)}{c^2 T^2} = \frac{y(n+1,m) - 2y(n,m) + y(n-1,m)}{T^2} \\
y(n, m+1) = y(n+1,m) + y(n-1,m) - y(n,m-1)
$$
We can solve this iteratively using some initial and boundary conditions, such as $y(t,0) = 0, y(0,x) = $ triangle function etc.