The easiest way to transform a state space model:
$$\dot{x} = Ax+Bu\\y = Cx + Du$$
Into a discrete state space equation $H(z)$ is to transform the state space equation and turn it to a discrete state space model:
$$\dot{x}_{k+1} = A_dx_k+B_du_k\\y = C_dx_k + D_du_k$$
Then compute the zeros, poles and gains, z, p, k. Use them to create the discrete transfer function:
$$H(z) = k\frac{\prod{z-z_i}}{\prod{z-p_i}}$$
Or in other words, compute the numerators and denominators.
So I have been using Octave, which is similar to MATLAB. Octave has a function named
expm() % Exponential matrix command
If I have a system matrix $A$ and a sample time $h$. To compute the discrete matrix $A_d$ I use this command:
$$A_d = \operatorname{expm}(Ah)$$
Very easy and fast. So I have now completed to turn $A$ into a discrete matrix. Now it's only $B$ matrix left. The definition of system matrix and signal matrix in a discrete state space model are:
$$A_d = e^{Ah}\\B_d = \int_0^h e^{A(h-t)}Bdt$$
The question is: How can I find the $B_d$ matrix? Finding the exponential matrix requries symbolic tools because the state transition matrix $e^{At}$ is not easy. There are servral ways. Here are the most used:
$$e^{At} = I + At + \frac{A^2t^2}{2!} + \frac{A^3t^3}{3!} \dots $$ $$e^{At} = L^{-1}(sI-A)^{-1}$$
The integral can be evaluated, but would then require you to also multiply by $A^{-1}$, which might not be possible if $A$ is singular. Instead a faster way of computing it can be used
$$ \begin{bmatrix} A_d & B_d \\ 0 & I \end{bmatrix} = e^{ \begin{bmatrix} A & B \\ 0 & 0 \end{bmatrix} h}, $$
with the size of the zeros such that the matrix in the exponent is square. Also use expm for the exponent.