diff --git a/Project.toml b/Project.toml
index c5b42e0..f08ee67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,7 @@ version = "0.2.1"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -14,7 +14,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 CSV = "~0.10, ~0.9, ~0.8, ~0.7"
-DifferentialEquations = "~7, ~6"
+OrdinaryDiffEq = "~6, ~5"
 JSON = "~0.21, ~0.20, ~0.19, ~0.18"
 julia = "^1"
 
diff --git a/src/QuantumAnnealing.jl b/src/QuantumAnnealing.jl
index cc79a1b..24a1342 100644
--- a/src/QuantumAnnealing.jl
+++ b/src/QuantumAnnealing.jl
@@ -4,7 +4,7 @@ module QuantumAnnealing
     import SparseArrays
     import JSON
     import Printf
-    import DifferentialEquations
+    import OrdinaryDiffEq
 
     include("base.jl")
     include("ising.jl")
diff --git a/src/simulate_de.jl b/src/simulate_de.jl
index 60a715b..5694964 100644
--- a/src/simulate_de.jl
+++ b/src/simulate_de.jl
@@ -1,6 +1,6 @@
 
 """
-simulator that uses the DifferentialEquations.jl differential equation solver
+simulator that uses the OrdinaryDiffEq.jl differential equation solver
 to solve the Schrodinger Equation, rather than the Magnus Expansion.  This
 simulator can run into issues at higher anneal times.  The API is the same as
 for the simulate function.
@@ -40,8 +40,8 @@ function simulate_de(ising_model, annealing_time, annealing_schedule, reltol; ab
     schrod_eq(state, time, s) = -im * time * H(s) * state
 
     s_range = (0.0, 1.0)
-    prob = DifferentialEquations.ODEProblem(schrod_eq, initial_state, s_range, annealing_time)
-    sol = DifferentialEquations.solve(prob; abstol=abstol, reltol=reltol, alg_hints=[:nonstiff], saveat=[1.0], kwargs...)
+    prob = OrdinaryDiffEq.ODEProblem(schrod_eq, initial_state, s_range, annealing_time)
+    sol = OrdinaryDiffEq.solve(prob, OrdinaryDiffEq.Tsit5(); abstol=abstol, reltol=reltol, saveat=[1.0], kwargs...)
 
 
     state = sol(1)
diff --git a/test/simulate.jl b/test/simulate.jl
index dcf8835..f6ab66a 100644
--- a/test/simulate.jl
+++ b/test/simulate.jl
@@ -459,40 +459,32 @@ end
     @testset "1 qubit, function schedule, constant terms" begin
         ρ = simulate_magnus_generic(one_spin_model, 1.0, AS_CIRCULAR, 100, 4, constant_field_x=[1], constant_field_z=[1])
         ρ_de = simulate_de(one_spin_model, 1.0, AS_CIRCULAR, 1e-6, constant_field_x=[1], constant_field_z=[1])
-
         @test isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "1 qubit, csv schedule, analytical solution" begin
         ρ = simulate_de(one_spin_model, 1.0, AS_CIRCULAR_pwl_csv_1000, 1e-6)
-
         @test isapprox(one_spin_ρ(1.0), ρ, atol=1e-6)
-        @test !isapprox(one_spin_ρ(1.0), ρ, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 1.0, AS_CIRCULAR, 100, 4)
         ρ_de = simulate_de(ising_model, 1.0, AS_CIRCULAR, 1e-6)
-
         @test isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule, long annealing time" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 1000.0, AS_CIRCULAR, 4000, 4)
-        ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-7)
-
+        ρ_de = simulate_de(ising_model, 1000.0, AS_CIRCULAR, 1e-8)
         @test isapprox(ρ, ρ_de, atol=1e-6)
-        @test !isapprox(ρ, ρ_de, atol=1e-7)
     end
 
     @testset "2 qubit, function schedule, adaptive tolerance" begin
         ising_model = Dict((1,) => -0.1, (2,) => -1, (1,2) => -1)
         ρ = simulate_magnus_optimized(ising_model, 100.0, AS_CIRCULAR, 500, 4)
         ρ_de = simulate_de(ising_model, 100.0, AS_CIRCULAR, silence=true)
-
-        @test isapprox(ρ, ρ_de, atol=1e-7)
-        @test !isapprox(ρ, ρ_de, atol=1e-8)
+        @test isapprox(ρ, ρ_de, atol=1e-6)
     end
 end