Hello Bernardo,
The different results here are maybe because of round-off errors between different machines or operating system, which are to be expected and which MATLAB cannot account for, as it is based on the machine itself. Because of the formulation of the ODE these numerical errors propagate and lead to the large differences seen in the final result. A simple view of these sort of problems can be seen by ODEs that require non-negativity:
ODE solvers perform numerical approximations of the derivatives and use those to update the solutions for small enough timesteps to get accurate solutions. This means that the errors locally can be small but when aggregated together over the entire time vector large differences can appear from expected solutions. For stiff problems in particular choosing the right set of solvers and constraints to get the ‘right’ answer is not obvious and can require some experimentation, as it may in your case. I recommend experimenting further with the formulation of your problem if reducing the error tolerance is not an effective workaround.
