Numerical problem with set of ODEs

Question

I have a set of ODEs:

$$\dot{x}_{i} = f_i(x_1, \ldots, x_N), ~ i \in \{1, 2, \ldots, N\}$$

This set of ODE has the following properties:

1. $\displaystyle\sum_{i=1}^N f_i = 0 ~ \forall x_1, \ldots, x_N$
2. If $\displaystyle x_i(0) \geq 0 ~ \forall i$ and $\displaystyle\sum_{i=1}^N x_i(0) = 1$ then $\displaystyle x_i(t) \geq 0 ~ \forall i, t > 0$ and $\displaystyle \sum_{i=1}^N x_i(t) = 1 ~\forall t > 0$
3. When $x_i(t) = 0$ or $x_i(t) = 1$ then $\dot{x}_i(t) = 0 ~ \forall x_j(t), i\neq j$.

Having said that, I have numerical problem while simulating these set of ODEs. In particular, I sometime get that for some $t>0$, I found $x_i(t) > 1$ or $x_i(t) < 0$. These errors can be small ($x_i(t) = 1 + 10^{-7}$ or $x_i(t) = 0 - 10^{-7}$) or even ''great'' ($x_i(t) = 1 + 10^{-3}$ or $x_i(t) = 0 - 10^{-3}$).

For simulation, I use standard MATLAB functions ode45, ode23 & co.

How can I prevent this behavior? If needed, I can add other details!