There is no advantage to using a better ODE solver on this example because the sensitivity to errors is so high. You can measure this using the uncertainty quantification methods in the Julia ODE solver suite:

https://docs.juliadiffeq.org/latest/analysis/uncertainty_qua...

where divergence on the Lorenz attactor tends to occur by t=80 or so even with accuracy of 1e-16.

But the funny thing about chaotic problems is that the shadowing theorem holds, which is:

>Although a numerically computed chaotic trajectory diverges exponentially from the true trajectory with the same initial coordinates, there exists an errorless trajectory with a slightly different initial condition that stays near ("shadows") the numerically computed one.

So you might as well just use Euler's method with high error, because it's backwards stable for this calculation, i.e. it gives a trajectory on the attractor, just the wrong trajectory. But since every method gives an O(1) error wrong trajectory after a short finite time, you might as well use the cheapest most error prone but convergent method.