Our package QuantumAnnealing solves ODEs that arise in design of quantum annealing algorithms. The package includes two solution methods `simulate`

(a hand crafted solution) and `simulate_de`

(a generic solution approach using DifferentialEquations).

Source code here.

We have noticed that there can be correctness issues with our `simulate_de`

implementation when running long time horizon problems, which makes me suspect we need to tweak the `solve`

settings we are using. To that end we are seeking advise on the following points,

- Are there changes we can make to improve the solution accuracy for long time horizons?
- Are there changes we can make to improve runtime performance?
- Are there any ways to avoid the precompile time of DifferentialEquations when running in CI? Compilation of DifferentialEquations is currently about 80% of our CI runtime.

Here is a MWEs to highlight the correctness issue. Note that `sum(z_measure_probabilities(ρ)`

should always be 1.0 in a perfect simulation.

```
using QuantumAnnealing
ising_model = Dict((1,) => 1.0, (1,2) => -1.0)
ρ = simulate_de(ising_model, 1.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 7.214235431263205e-9
ρ = simulate_de(ising_model, 10.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 8.508320654354584e-6
ρ = simulate_de(ising_model, 100.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.0003582731716551546
ρ = simulate_de(ising_model, 1000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.003550520566418469
ρ = simulate_de(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.03483614681737934
ρ = simulate_de(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.2990146971495551
ρ = simulate(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.3211653993039363e-14
ρ = simulate(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.0436096431476471e-13
ρ_target = simulate(ising_model, 10000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 10000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.0348400079153087
ρ_target = simulate(ising_model, 100000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 100000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.2990175297085176
```

CC @zmorrell