It's not a complete answer but I hope that it could be of some help.
You could rewrite the first system as
$$
\begin{cases}
P(n) = K_P E(n) \\
I(n) = I(n-1) + \frac{K_P}{T_I} E(n) \Delta t \\
D(n) = K_P T_D \frac{E(n) - E(n-1)}{\Delta t}
\end{cases}
$$
Where \$E(n) = G(n) - target(n)\$ and \$\Delta t\$ is your sampling interval. Note that \$T_D\$ and \$T_I\$ are not defined as gains. \$K_I = \frac{K_P}{T_I}\$ and \$K_D = K_P T_I\$ are respectively the integral gain and the derivative gain.
Now you can rewrite the system as a single function of the error.
$$
PID(n) = P(n) + I(n) + D(n)
$$
$$
I(n-1) = PID(n-1) - P(n-1) - D(n-1) \\
= PID(n-1) - K_P E(n-1) - K_P T_D \frac{E(n-1) - E(n-2)}{\Delta t}
$$
$$
PID(n) = K_P E(n) + PID(n-1) - K_P E(n-1) - K_P T_D \frac{E(n-1) - E(n-2)}{\Delta t} + \frac{K_P}{T_I} E(n) \Delta t + K_P T_D \frac{E(n) - E(n-1)}{\Delta t} \\
= PID(n-1) + K_P \left(\left(1 + \frac{\Delta t}{T_I} + \frac{T_D}{\Delta t} \right)E(n) - \left(1 + 2\frac{T_D}{\Delta t} \right)E(n-1) + \frac{T_D}{\Delta t} E(n-2) \right)
$$
The second one is a bit more complex to rewrite as a single equation but you can do it in a similar way. The result should be
$$
R(n) = K_1 R(n-1) - (\gamma K_0 + K_2) R(n-2) + (1+\gamma) (PID(n) - K_1 PID(n-1) + K_2 PID(n-2))
$$
Now you only need to substitute the equation of the PID in order to obtain the equation of the regulator as function of the error.
On the integration mentioned in the pdf
To check for the zero supposedly sitting at \$-2+0i\$ , the integration required is
\$\int_0^\infty k e^{-1t}sin(3t+\phi) e^{-(-2+0 i) t} dt\$.
But that value of \$s = -2+0 i\$ is outside the Region of convergence (ROC) for the given signal. You already saw that it is not converging. The Laplace transform exists only to the right side of the right most pole which in this case is \$-1\pm 3i\$.
See slide #10 of this presentation and
Wikipedia
In the pdf that you have linked, the zero happens to be in the ROC. It is to the right of the right most pole. Figure is shown below.
You can adjust the spring mass damper example in your case so that the zero is to the right of the right most pole and try again to see if you are getting the result you expect (convergence of the integral).
Also, I still think there may be some confusion regarding the zeros of a specific signal (solved using initial conditions) and that of a system. The pdf does not mention initial conditions; they would have assumed IC as zero. However, non-zero initial conditions appear in your spring mass example.
On the concept of Zeroes for a spring mass system
A spring mass system with one spring and mass can be written in many ways (Your spring mass system input has been defined as zero, making the determination of zeroes of the system difficult).
Each way changes the zeroes of the system. Some examples are given below.
system (A)
Input : force f(t)
Output : displacement of mass y(t)
Net force on the mass : f(t) + spring force
Zeros of system : none (or two zeroes at \$s = \infty i\$)
Intuition : When force is varying with infinite frequency, the mass doesn't move due to its inertia.
Expressions
\$
\frac{d^2 y(t)}{dt^2} = (f(t) - Ky(t))/M
\$
\$
\frac{Y(s)}{F(s)} = \frac{1/M}{s^2 + K/M}
\$
system(B)
Input : force f(t)
Output : a force measurement sensor output m(t) (which happens to depend on the displacement and velocity of the moving mass)
Net force on the mass : f(t) + spring force (\$K_1\$ can be combined with \$K\$) + dash pot force
Zeros of system : one zero at \$s = -K_1/C_1 + 0i\$
Intuition : When \$y(t) > 0\$ spring \$K_1\$ excerts a right-ward force on the sensor.
At the same time if \$dy(t)/dt < 0\$, then the dash pot excerts a left-ward force on the sensor.
If \$K_1y(t) = -C_1 dy(t)/dt\$, the sensor senses zero force; i.e output of system is zero even when the internal states are not zero.
i.e. a zero of the system.
The condition \$K_1y(t) = -C_1 dy(t)/dt\$ means that the displacement is expressed as
\$y(t) = e^{-t K_1/C_1}\$; Thus connecting the whole thing to an exponential signal.
Expressions
\$
\frac{d^2 y(t)}{dt^2} = (f(t) - (K+K_1) y(t))/M - C_1 \frac{d y(t)}{dt}
\$
\$
\frac{Y(s)}{F(s)} = \frac{1/M}{s^2 + C_1 s + (K+K_1)/M}
\$
\$
m(t) = K_1 y(t) + C_1 \frac{d y(t)}{dt}
\$
\$
\frac{M(s)}{F(s)} = \frac{M(s)}{Y(s)} \frac{Y(s)}{F(s)} =
\frac{(s C_1 + K_1) 1/M}{s^2 + C_1 s + (K+K_1)/M}
\$
system (C)
Input : displacement x(t)
Output : displacement of mass y(t)
Net force on the mass : spring force (spring extension is \$x(t)-y(t)\$)
Zeros of system : none (or two zeroes at \$s = \infty i\$)
Expressions
\$
\frac{d^2 y(t)}{dt^2} = (x(t)-y(t))K/M
\$
\$
\frac{Y(s)}{X(s)} = \frac{K/M}{s^2 +K/M}
\$
system (D)
Input : displacement x(t)
Output : displacement of mass y(t)
Net force on the mass : spring force + dash pot force
Zeros of system : one zero at \$s = -K/C + 0i\$
Intuition : When \$y(t) = 0, Kx(t)= -Cdx(t)/dt\$, the net force on the mass is zero. Hence the mass remains stationary;
i.e. output of system is zero, i.e. a zero of the system.
Expressions
\$
\frac{d^2 y(t)}{dt^2} = (x(t)-y(t))K/M + \frac{d(x(t)-y(t))}{dt}C/M
\$
\$
\frac{Y(s)}{X(s)} = \frac{s C/M + K/M}{s^2 + s C/M + K/M}
\$
It can be noted that in all cases, the initial condition of the system is not mentioned / assumed zero. The zeroes of the system depends strongly on the way in which input and output of the system are defined as well as the dynamics of the system.
Best Answer
Consider the following 2nd order system:
$$a\ddot y(t) + b\dot y(t) + cy(t) = x(t)$$
We can write this as
$$\mathbf Dy(t) = x(t) $$
where \$\mathbf D\$ is the 2nd order operator
$$\mathbf D = a\frac{d^2}{dt^2} + b\frac{d}{dt} + c$$
We know that \$e^{st}\$ is an eigenfunction of the operator \$\mathbf D\$ with eigenvalue \$\lambda = as^2 + bs + c\$
$$\mathbf De^{st} = \lambda e^{st} = (as^2 + bs + c)e^{st}$$
thus there are an infinity of independent eigenfunctions for this linear time invariant differential operator.
Next, consider the equation
$$\mathbf Dy(t) = 0$$
Now we have a constraint - we seek the specific eigenfunction(s) with eigenvalue \$\lambda =0\$.
As is well known, we proceed as follows:
$$\lambda = as^2 + bs + c = 0$$
which has two independent solutions
$$s = s_{0\pm} = -\frac{b}{2a} \pm \sqrt{\frac{b^2}{(2a)^2} - \frac{c}{a}}$$
These particular eigenfunctions are clearly special in that they are the modes of the system - the form of the system output when there is no input.