Given the differential equation with variables $x_{1}(t)$, $x_2(t)$, $u_{1}(t)$:
$$\frac{dx_{1}}{dt}=\sqrt{k\cdot{x_1}}+\sqrt{k_{1}\cdot x_{2}}+{u_1}$$ I have linearised it around the equilibrium point ($\frac{dh_{1}(t)}{dt}=0$) via a first order Taylor expansion to get
$$\frac{dx_{1}}{dt}=\frac{k}{2\sqrt{k\cdot x_{1\hspace{1mm}equilibria}}}\cdot \Delta{x_1}+\frac{k_1}{2\sqrt{k_1\cdot x_{2\hspace{1mm}equilibria}}}\cdot \Delta x_{2}+\Delta u_{1}$$ Now, by making (just for simplicity) $$a=\frac{k}{2\sqrt{k\cdot x_{1\hspace{1mm}equilibria}}}, b=\frac{k_1}{2\sqrt{k_1\cdot x_{2\hspace{1mm}equilibria}}}$$
the equation becomes
$$\frac{dx_{1}}{dt}=a\cdot \Delta{x_1}+b\cdot \Delta x_{2}+\Delta u_{1}$$
Until this point I had no trouble, but I have to discretise the linearised differential equation, in order to apply the $Z$-Transform, with sampling time $T=1s$.
My thoughts are that $\frac{dx_{1}}{dt}$ can be approximated with finite differences, thus resulting in
$$\frac{dx_{1}}{dt}=\frac{x_1((n+1)T)-x_1(nT)}{T}=\Delta x_1(n) =a\cdot \Delta{x_1}+b\cdot \Delta x_{2}+\Delta u_{1} $$
However, I am not sure on how to proceed now when it comes to discretising the other deltas. If it helps, it seems that the expected result is
$$\Delta x_1(n+1) =a\cdot \Delta{x_1}(n)+b\cdot \Delta x_{2}(n)+\Delta u_{1}(n)$$
Any ideas or thoughts on how can the linearised equation be discretised?
Normally if you would have some non-linear time-invariant multidimensional differential equation with input vector $\vec{u}(t)$ of the form,
$$ \frac{d}{dt}\vec{x}(t) = \vec{f}(\vec{x}(t), \vec{u}(t)), \tag{1} $$
with,
$$ \vec{x}(t) = \begin{bmatrix} x_1(t) & x_2(t) & \cdots & x_n(t) \end{bmatrix}^T, \tag{2} $$
$$ \vec{u}(t) = \begin{bmatrix} u_1(t) & u_2(t) & \cdots & u_m(t) \end{bmatrix}^T, \tag{3} $$
$$ \vec{f}(\vec{x}, \vec{u}) = \begin{bmatrix} f_1(\vec{x}, \vec{u}) & f_2(\vec{x}, \vec{u}) & \cdots & f_n(\vec{x}, \vec{u}) \end{bmatrix}^T, \tag{4} $$
you could linearize it around an equilibrium point $\vec{x}_{eq}$ (such that $\vec{f}(\vec{x}_{eq}, \vec{0})=\vec{0}$). Such a linearization would yield a differential equation of the following form,
$$ \vec{z}(t) = \vec{x}(t) - \vec{x}_{eq}, \tag{5} $$
$$ \frac{d}{dt}\vec{z}(t) = A\,\vec{z}(t) + B\,u(t), \tag{6} $$
where $A$ is a $n$ by $n$ matrix and $B$ a $n$ by $m$ matrix. Matrix $A$ will be equal to the Jacobian evaluated at the linearization point, which is equivalent to your first order Taylor expansion, which can be found with,
$$ A = \left.\begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \cdots & \frac{\partial f_n}{\partial x_n} \end{bmatrix}\right|_{(\vec{x},\vec{u})=(\vec{x}_{eq},\vec{0})}. \tag{7} $$
Matrix $B$ can be found in a similar way,
$$ B = \left.\begin{bmatrix} \frac{\partial f_1}{\partial u_1} & \cdots & \frac{\partial f_1}{\partial u_m} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial u_1} & \cdots & \frac{\partial f_n}{\partial u_m} \end{bmatrix}\right|_{(\vec{x},\vec{u})=(\vec{x}_{eq},\vec{0})}. \tag{8} $$
The solution for $\vec{z}(t)$ subjected to equation $(6)$ and initial condition $\vec{z}(t_0)=\vec{z}_0$ for $t\geq t_0$ can be found to be,
$$ \vec{z}(t) = e^{A\,(t-t_0)}\,\vec{z}_0 + \int_{t_0}^t e^{A\,(t - \tau)}\,B\,\vec{u}(\tau)\,d\tau. \tag{9} $$
Now if you would like to simulate this system in steps of $T$ and assume that,
$$ \vec{u}(t) = \vec{u}(k\,T), \quad k\,T \leq t < (k+1)\,T, \quad \forall\ \ k\in\mathbb{Z}, \tag{10} $$
then for $t_0=k\,T$ and $t=(k+1)\,T$, equation $(9)$ can be written as,
$$ \vec{z}((k+1)\,T) = \underbrace{e^{A\,T}}_{A_d}\,\vec{z}(k\,T) + \underbrace{\int_{k\,T}^{(k+1)\,T}\mkern-18mu e^{A\,((k+1)\,T - \tau)}\,d\tau\,B}_{B_d = A^{-1}\left(A_d-I\right)B}\,\vec{u}(k\,T). \tag{11} $$
This is equivalent to the discrete representation,
$$ \vec{z}[k+1] = A_d\,\vec{z}[k] + B_d\,\vec{u}[k]. \tag{12} $$
For your system you only gave $f_1(\vec{x},\vec{u})$, so you would also need $f_2(\vec{x},\vec{u})$ in order to complete the total linearization and subsequently the discretization.
You could also use your approach, but it will be less accurate. Namely this would only do a first order approximation of the matrix exponentiation and will only yield approximately the same results for very small $T$,
$$ e^{A\,T} = I + A\,T + \frac12 A^2\,T^2 + \frac{1}{3!} A^3\,T^3 + \cdots. \tag{13} $$