I shall try to get to the point, so let me know if there is something left and you need more details.
I am solving a couple of equations that are not coupled explicitly, but their corresponding unknown variables, say $x$ and $y$ must satisfy a differential equation:
$\dot x = x + y,$
where dots denote derivative with respect to an independent variable, say $t$.
The equation for $x$ is of second-order, so one gets $x$ and $\dot x$ from it, and one can check if the equation above is consistently satisfied. However (see the plot attached), it turns out that, no matter which integration method I use from SciPy (the ones that are already implemented), the equality above stops being satisfied at some point. This is due to the fact that $x$ and $y$ cancel each other up to a very high precision, which seems not to be attained by none of the methods provided by SciPy (I have checked this by taking every method and by lowering the absolute and relative tolerance as much as possible. In the plot attached, the method employed is DOP853, which is supposed to be very useful when very low tolerances are required).
My question is if you know of any way to improve the accuracy so the cancelation gets more precise (I would like that the equation will be conveniently satisfied all the way through the entire calculation). The only parameters I have changed so far were the relative and absolute tolerances (and of course the different methods at our disposal). Is there any parameter that I am missing and that might be useful for that?