Problem in verifying zero blocking property of a zero of a multiple input multiple output system
The transmission zeros of a MIMO system are obtained by solving

\begin{align}
\dot{x}^* = A\,x^* + B\,u^*, \\
y = C\,x^* + D\,u^*,
\end{align}

for $x^*$ and $u^*$ such that $y = 0$. By using the Laplace variable $s$ for the time derivative that system of equation can also be written as

$$
\begin{bmatrix}
0 \\ 0
\end{bmatrix} =
\begin{bmatrix}
A - s\,I &amp; B \\ C &amp; D
\end{bmatrix}
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

which can be seen as the generalized eigenvalue problem $\mathcal{A}\,v = \lambda\,\mathcal{B}\,v$, with $\lambda$ equivalent to $s$,

$$
\mathcal{A} = 
\begin{bmatrix}
A &amp; B \\ C &amp; D
\end{bmatrix}, \quad
\mathcal{B} = 
\begin{bmatrix}
I &amp; 0 \\ 0 &amp; 0
\end{bmatrix}, \quad
v = 
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix}.
$$

I believe that tzero in matlab also solves such generalized eigenvalue problem to obtain the transmission zeros (however by using eig yourself the associated eigenvectors are more accurate than the nullspace calculations). The associated solution in time for both $x(t)$ and $u(t)$ can be obtained using

$$
\begin{bmatrix}
x(t) \\ u(t)
\end{bmatrix} = e^{\lambda\,t}\,
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

so you can't just pick any value for $\lambda$ when generating such solution. You could now simulate the system using $x(0)=x^*$ as initial condition and the expression for $u(t) = e^{\lambda\,t}u^*$ with lsim. However, I believe that slim uses zero-order-hold for the input so the actual applied input will deviate slightly causing the output $y$ to deviate away from zero using that approach. Instead, if you look at the analytical solution, both the output and the error in state dynamics should remain practically to zero (for non-positive $\lambda$):


  n = length(A);
  m = length(D);
  [V,L] = eig([A B; C D], diag([ones(1,n) zeros(1,m)]));
  ind = 1;    % Index of associated eigenvalue (avoid infinite eigenvalues)
  lambda = L(ind,ind);
  x0 = V(1:n, ind);
  u0 = V(1+n:n+m, ind);  
  
  N = 1e3;
  t = linspace(0, 5, N);
  x = @(t) x0 * exp(lambda* t);
  u = @(t) u0 * exp(lambda* t);
  y = @(t) C * x(t) + D * u(t);
  e = @(t) A * x(t) + B * u(t) - lambda* x(t);  
  
  Y = zeros(n, N);
  E = zeros(m, N);
  for k = 1 : N
      Y(:,k) = y(t(k));
      E(:,k) = e(t(k));
  end
  figure, plot(t, E)
  figure, plot(t, Y)