Due to round-off errors and local discretization error the numerical and mathematical behavior is expected to be different. A simply example is the differential quotient as approximation of the derivative:
If h is large, the linearization error is large. If h is small, the difference of the function evaluations can be such tiny, that the cancellation error dominates the result. Therefore the numerically determined derivative has a limited accuracy only, as a rule of thumb in the magnitude of sqrt(eps). Equivalent arguments appear for the two-sided differential quotient and any other scheme also.
Because all numerical ODE solvers are based on evaluating linearized derivatives, the limited accuracy is a fundamental effect, which cannot be avoided. Reducing the relative and absolute tolerances helps, but only until a certain limit, when the cancellation error dominates the linearization error and the accumulation of the round-off errors gets more important due to the higher number of steps.
However, an unexpected instability can be caused by a modelling error also. Did you validate the model exhaustively? How did you check the dynamic instability? Is the model continuously differentiable (I ask this, because in the past I saw so many users in this forum, who used integrands with discontinuities)?