Is there any alternative way to simulate a continuous time model instead of discrete?

Question

Normally I use the exponential matrix method to turn continuous time model to discrete time model. In this case, it's SS-models:

$$\dot x = Ax + Bu + Eu \\ y = Cx + Du$$

$$\begin{bmatrix} A_d & B_d &E_d \\ 0 &I &I \end{bmatrix} = e^{\begin{bmatrix} A & B & E\\ 0 & 0 & 0 \end{bmatrix}}T_s$$

Where $T_s$ is the sampling time.

But I have a problem with it comes to this and it's because of the simulation of nonlinear systems.

If I have a nonlinear model:

$$\dot x = f(x, u)$$

And I want to linearize it in the state $x_0 = 0$. I know that it will be nothing left of the dynamics. Assume that we linearize at the initial state and we have some damping in our system, then the damping factor is going to be zero.

I going to do nonlinear numerical simulation of a discrete system. So my procedure is going to be:

1. Linearize the system in state $x$. E.g If we have the damping $Bx_2^2$ the it would be linearized to $2Bx_2$.
2. Turn it to discrete
3. Do the simulation and find state $x$
4. Jump to stage 1.

A lot of times, my simulation breaks. To do a simulation, I need to linearize the system first and then turn it to a discrete form, as I have been learned.

Question:

Is there any simpler way to simulate a nonlinear continuous time model? Is there an alternative way to turn the model to discrete?

Nonlinear simulation seems to be so sensitive to against everything.

EDIT:

Here is the code. Assume that you know A, B and sample time h.

% Compute sizes
a1 = size(A,2) + size(B,2) - size(A,1);
b1 = size(A,2);
a2 = size(A,2) + size(B,2) - size(B,1);
b2 = size(B,2);
% Compute square matrix
M = [A B; zeros(a1, b1)  zeros(a2, b2)];
M = expm(M*h);
% Find the discrete matrecies
Ad = M(1:size(A,1), 1:size(A,2));
Bd = M(1:size(B,1), (size(A,2) + 1):(size(A,2) + size(B,2)));