diff --git a/test/weak_convergence/PL1WM.jl b/test/weak_convergence/PL1WM.jl
index 11522f0d6..315b8920c 100644
--- a/test/weak_convergence/PL1WM.jl
+++ b/test/weak_convergence/PL1WM.jl
@@ -7,18 +7,9 @@ import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
 end
@@ -50,15 +41,14 @@ ensemble_prob = EnsembleProblem(prob;
         output_func = (sol,i) -> (h1(sol[end]),false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
 
-#using Plots; convergence_plot = plot(dts, errors, xaxis=:log, yaxis=:log)
-#savefig(convergence_plot, "PL1WM-scalar.png")
-println("PL1WM:", m)
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.34 # order is 2.34
+println("PL1WM:", sim.𝒪est[:weak_final])
 
 
 """
@@ -87,14 +77,13 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("PL1WM:", m)
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.34 # order is 2.34
+println("PL1WM:", sim.𝒪est[:weak_final])
 
 """
  Test non-commutative noise SDEs (iip)
@@ -144,15 +133,13 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-u₀[1]*exp(2*1.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.33
-
-println("PL1WM:", m)
-
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀[1]*exp(2*1.0)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.33 # order is 1.6748033428458136
+println("PL1WM:", sim.𝒪est[:weak_final])
 
 
 """
@@ -186,14 +173,13 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("PL1WM:", m)
-
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("PL1WM:", sim.𝒪est[:weak_final])
 
 
 """
@@ -217,15 +203,26 @@ ensemble_prob = EnsembleProblem(prob;
         output_func = (sol,i) -> (sol[end],false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-_solutions1 = @time generate_weak_solutions(ensemble_prob, PL1WMA(), dts, numtraj, ensemblealg=EnsembleThreads())
 
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-@test minimum(_solutions .≈ _solutions1)
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+    )
+
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("PL1WM:", sim.𝒪est[:weak_final])
 
-println("PL1WM:", m)
+sim1 = test_convergence(dts,ensemble_prob,PL1WMA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+    )
+
+@test abs(sim1.𝒪est[:weak_final]-2) < 0.3
+println("PL1WMA:", sim1.𝒪est[:weak_final])
+
+@test minimum(sim.solutions .≈ sim1.solutions)
 
 #inplace
 
@@ -236,13 +233,27 @@ gadd!(du,u,p,t) = @.(du = p[2])
 
 prob = SDEProblem(fadd!,gadd!,u₀,tspan,p)
 ensemble_prob = EnsembleProblem(prob;
-        output_func = (sol,i) -> (sol[end],false),
+        output_func = (sol,i) -> (sol[end][1],false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, PL1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-_solutions1 = @time generate_weak_solutions(ensemble_prob, PL1WMA(), dts, numtraj, ensemblealg=EnsembleThreads())
 
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-@test minimum(_solutions .≈ _solutions1)
+
+sim = test_convergence(dts,ensemble_prob,PL1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("PL1WM:", sim.𝒪est[:weak_final])
+
+sim1 = test_convergence(dts,ensemble_prob,PL1WMA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+    )
+
+@test abs(sim1.𝒪est[:weak_final]-2) < 0.3
+println("PL1WMA:", sim1.𝒪est[:weak_final])
+
+@test minimum(sim.solutions .≈ sim1.solutions)
diff --git a/test/weak_convergence/SIE_SME.jl b/test/weak_convergence/SIE_SME.jl
index 94c2d07af..88afc2e78 100644
--- a/test/weak_convergence/SIE_SME.jl
+++ b/test/weak_convergence/SIE_SME.jl
@@ -7,18 +7,9 @@ import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
 end
@@ -50,42 +41,38 @@ ensemble_prob = EnsembleProblem(prob;
         output_func = (sol,i) -> (h1(sol[end]),false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-#using Plots; convergence_plot = plot(dts, errors, xaxis=:log, yaxis=:log)
-#savefig(convergence_plot, " SIEA.png")
-println("SIEA:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEA:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
 
-println("SIEB:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEB:", m)
+sim = test_convergence(dts,ensemble_prob,SIEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
+println("SIEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
+println("SMEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SIEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
+println("SIEB:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
+println("SMEB:", sim.𝒪est[:weak_final])
 
 """
  Test Scalar SDEs (iip)
@@ -113,41 +100,38 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
+sim = test_convergence(dts,ensemble_prob,SIEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
+println("SIEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
+println("SMEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SIEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
+println("SIEB:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀.*exp(1.0*(p[1]))
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
+println("SMEB:", sim.𝒪est[:weak_final])
 
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SIEA:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEA:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SIEB:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors =  [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEB:", m)
 
 """
  Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
@@ -180,35 +164,34 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SIEA:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEA(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEA:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SIEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SIEB:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, SMEB(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("SMEB:", m)
+sim = test_convergence(dts,ensemble_prob,SIEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("SIEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEA(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("SMEA:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SIEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("SIEB:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,SMEB(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("SMEB:", sim.𝒪est[:weak_final])
diff --git a/test/weak_convergence/srk_weak_final.jl b/test/weak_convergence/srk_weak_final.jl
index 51089f46a..de4b3a7f6 100644
--- a/test/weak_convergence/srk_weak_final.jl
+++ b/test/weak_convergence/srk_weak_final.jl
@@ -3,29 +3,14 @@
  DRI1, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM
 """
 
-
 import Statistics # for mean values of trajectories
 import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-#using Plots
-
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],adaptive=false,
-      save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,
-      trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
-
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
 end
@@ -57,116 +42,111 @@ ensemble_prob = EnsembleProblem(prob;
         output_func = (sol,i) -> (h1(asinh(sol[end])),false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1(), dts, numtraj, ensemblealg=EnsembleThreads())
 
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
 
-#convergence_plot = plot(dts, errors, xaxis=:log, yaxis=:log)
-#savefig(convergence_plot, "DRI1-CPU-final-"*string(numtraj)*".pdf")
-println("DRI1:", m)
+sim = test_convergence(dts,ensemble_prob,DRI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("DRI1:", sim.𝒪est[:weak_final])
+
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI1:", m)
+sim = test_convergence(dts,ensemble_prob,RI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI1:", sim.𝒪est[:weak_final])
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI3(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI3:", m)
+sim = test_convergence(dts,ensemble_prob,RI3(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI3:", sim.𝒪est[:weak_final])
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI5(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI5:", m)
-
+sim = test_convergence(dts,ensemble_prob,RI5(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI5:", sim.𝒪est[:weak_final])
 
 numtraj = Int(7e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI6(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI6:", m)
+sim = test_convergence(dts,ensemble_prob,RI6(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI6:", sim.𝒪est[:weak_final])
 
 numtraj = Int(1e3)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-1) < 0.3
-
-println("RDI1WM:", m)
-
+sim = test_convergence(dts,ensemble_prob,RDI1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-1) < 0.3
+println("RDI1WM:", sim.𝒪est[:weak_final])
 
 numtraj = Int(7e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI2WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI2WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI3WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-numtraj = Int(3e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-println("RDI3WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI4WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
 
-println("RDI4WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI2WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI2WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI3WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI3WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI4WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.568 with 3e5 trajectories
+println("RDI4WM:", sim.𝒪est[:weak_final])
 
 """
  Test Scalar SDEs (iip)
@@ -195,120 +175,119 @@ Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
 
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
+sim = test_convergence(dts,ensemble_prob,DRI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("DRI1:", sim.𝒪est[:weak_final])
 
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1NM(), dts, numtraj, ensemblealg=EnsembleThreads())
 
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("DRI1NM:", m)
+sim = test_convergence(dts,ensemble_prob,DRI1NM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("DRI1NM:", sim.𝒪est[:weak_final])
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI1:", m)
+sim = test_convergence(dts,ensemble_prob,RI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI1:", sim.𝒪est[:weak_final])
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI3(), dts, numtraj, ensemblealg=EnsembleThreads())
+sim = test_convergence(dts,ensemble_prob,RI3(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI3:", sim.𝒪est[:weak_final])
 
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI3:", m)
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI5(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI5:", m)
-
+sim = test_convergence(dts,ensemble_prob,RI5(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI5:", sim.𝒪est[:weak_final])
 
 numtraj = Int(7e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI6(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI6:", m)
+sim = test_convergence(dts,ensemble_prob,RI6(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI6:", sim.𝒪est[:weak_final])
 
 numtraj = Int(1e3)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-1) < 0.3
+sim = test_convergence(dts,ensemble_prob,RDI1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-1) < 0.3
+println("RDI1WM:", sim.𝒪est[:weak_final])
 
-println("RDI1WM:", m)
 
 numtraj = Int(7e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI2WM(), dts, numtraj, ensemblealg=EnsembleThreads())
+sim = test_convergence(dts,ensemble_prob,RDI2WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI2WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI3WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI3WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI4WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=0.0
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI4WM:", sim.𝒪est[:weak_final])
 
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI2WM:", m)
-
-numtraj = Int(3e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI3WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI3WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI4WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI4WM:", m)
 
 """
  Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
@@ -341,48 +320,46 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("DRI1:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1NM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("DRI1NM:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI1:", m)
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RI3(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI3:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RI5(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI5:", m)
+sim = test_convergence(dts,ensemble_prob,DRI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
+println("DRI1:", sim.𝒪est[:weak_final])
+
+
+sim = test_convergence(dts,ensemble_prob,DRI1NM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
+println("DRI1NM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 3.05
+println("RI1:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RI3(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.77
+println("RI3:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RI5(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI5:", sim.𝒪est[:weak_final])
 
 
 numtraj = Int(1e5)
@@ -390,13 +367,13 @@ seed = 70
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RI6(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.45
-
-println("RI6:", m)
+sim = test_convergence(dts,ensemble_prob,RI6(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.45 # order is 1.55
+println("R6:", sim.𝒪est[:weak_final])
 
 
 numtraj = Int(1e3)
@@ -404,13 +381,13 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-1) < 0.45
-
-println("RDI1WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-1) < 0.45 # order is 1.44
+println("RDI1WM:", sim.𝒪est[:weak_final])
 
 
 numtraj = Int(1e5)
@@ -418,13 +395,13 @@ seed = 70
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI2WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.45
-
-println("RDI2WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI2WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.45 # order is 1.55
+println("RDI2WM:", sim.𝒪est[:weak_final])
 
 
 numtraj = Int(5e4)
@@ -432,18 +409,18 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI3WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI3WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI4WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI4WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI3WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.84
+println("RDI3WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI4WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order is 2.91
+println("RDI4WM:", sim.𝒪est[:weak_final])
diff --git a/test/weak_convergence/srk_weak_final_non_diagonal.jl b/test/weak_convergence/srk_weak_final_non_diagonal.jl
index e3897b642..9ad481a8a 100644
--- a/test/weak_convergence/srk_weak_final_non_diagonal.jl
+++ b/test/weak_convergence/srk_weak_final_non_diagonal.jl
@@ -9,22 +9,9 @@ import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-#using Plots
-
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],adaptive=false,
-      save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,
-      trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
 
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
@@ -65,61 +52,61 @@ Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
 
-_solutions = @time generate_weak_solutions(ensemble_prob, DRI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("DRI1:", m)
+sim = test_convergence(dts,ensemble_prob,DRI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("DRI1:", sim.𝒪est[:weak_final])
 
 numtraj = Int(2e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
-_solutions = @time generate_weak_solutions(ensemble_prob, RI1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI1:", m)
+sim = test_convergence(dts,ensemble_prob,RI1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order 3.43
+println("RI1:", sim.𝒪est[:weak_final])
 
 numtraj = Int(2e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
-_solutions = @time generate_weak_solutions(ensemble_prob, RI3(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI3:", m)
+sim = test_convergence(dts,ensemble_prob,RI3(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order 2.46
+println("RI3:", sim.𝒪est[:weak_final])
 
 numtraj = Int(2e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
-_solutions = @time generate_weak_solutions(ensemble_prob, RI5(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI5:", m)
+sim = test_convergence(dts,ensemble_prob,RI5(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order 2.57
+println("RI5:", sim.𝒪est[:weak_final])
 
 numtraj = Int(2e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
-_solutions = @time generate_weak_solutions(ensemble_prob, RI6(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RI6:", m)
+sim = test_convergence(dts,ensemble_prob,RI6(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RI6:", sim.𝒪est[:weak_final])
 
 
 numtraj = Int(1e3)
@@ -127,39 +114,39 @@ seed = 10
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI1WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-1) < 0.3
-
-println("RDI1WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI1WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI1WM:", sim.𝒪est[:weak_final])
 
 numtraj = Int(1e5)
 seed = 10
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI2WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI2WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI3WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI3WM:", m)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RDI4WM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-exp(-3.0)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test -(m-2) < 0.3
-
-println("RDI4WM:", m)
+sim = test_convergence(dts,ensemble_prob,RDI2WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test -(sim.𝒪est[:weak_final]-2) < 0.3 # order 2.517769274990593
+println("RDI2WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI3WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI3WM:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RDI4WM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=exp(-3.0)
+)
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RDI4WM:", sim.𝒪est[:weak_final])
diff --git a/test/weak_convergence/weak_strat.jl b/test/weak_convergence/weak_strat.jl
index 863bef00e..91f1486c0 100644
--- a/test/weak_convergence/weak_strat.jl
+++ b/test/weak_convergence/weak_strat.jl
@@ -5,24 +5,14 @@
  NON, COM, NON2
 """
 
-
 import Statistics # for mean values of trajectories
 import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],adaptive=false,save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
 end
@@ -54,75 +44,58 @@ ensemble_prob = EnsembleProblem(prob;
         output_func = (sol,i) -> (h1(sol[end]),false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) - u₀*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-println("RS1:", m)
-
-
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())
 
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) - u₀*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.37
 
-println("RS2:", m)
+sim = test_convergence(dts,ensemble_prob,RS1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RS1:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RS2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.36 #order is 2.36
+println("RS2:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,NON(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON:", sim.𝒪est[:weak_final])
 
+numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
+sim = test_convergence(dts,ensemble_prob,NON2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON2:", sim.𝒪est[:weak_final])
 
-prob = SDEProblem(f,g,u₀,tspan, p)
-ensemble_prob = EnsembleProblem(prob;
-        output_func = (sol,i) -> (h1(sol[end]),false),
-        prob_func = prob_func
-        )
-
-
-_solutions = @time generate_weak_solutions(ensemble_prob, NON(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON:", m)
 
 dts = 1 .//2 .^(4:-1:0)
-
 numtraj = Int(1e5)
 seed = 10
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, COM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("COM:", m)
-
-dts = 1 .//2 .^(5:-1:0)
-numtraj = Int(5e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, NON2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON2:", m)
+sim = test_convergence(dts,ensemble_prob,COM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("COM:", sim.𝒪est[:weak_final])
 
 
 """
@@ -138,56 +111,50 @@ numtraj = Int(1e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
+dts = 1 .//2 .^(5:-1:0)
 
 prob = SDEProblem(f!,g!,[u₀],tspan, p)
 ensemble_prob = EnsembleProblem(prob;
-        output_func = (sol,i) -> (h1(sol[end]),false),
-        prob_func = prob_func
-        )
-_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts,
-  numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-println("RS1:", m)
-
-
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-prob = SDEProblem(f!,g!,[u₀],tspan, p)
-ensemble_prob = EnsembleProblem(prob;
-        output_func = (sol,i) -> (h1(sol[end]),false),
+        output_func = (sol,i) -> (h1(sol[end][1]),false),
         prob_func = prob_func
         )
-_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts,
-  numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.37
 
-println("RS1:", m)
+sim = test_convergence(dts,ensemble_prob,RS1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RS1:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RS2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.36 #order is 2.36
+println("RS2:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,NON(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON:", sim.𝒪est[:weak_final])
 
-numtraj = Int(1e5)
+numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
+sim = test_convergence(dts,ensemble_prob,NON2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON2:", sim.𝒪est[:weak_final])
 
-prob = SDEProblem(f!,g!,[u₀],tspan, p)
-ensemble_prob = EnsembleProblem(prob;
-        output_func = (sol,i) -> (h1(sol[end]),false),
-        prob_func = prob_func
-        )
-_solutions = @time generate_weak_solutions(ensemble_prob, NON(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON:", m)
 
 dts = 1 .//2 .^(4:-1:0)
 numtraj = Int(1e5)
@@ -195,28 +162,13 @@ seed = 10
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, COM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("COM:", m)
-
-dts = 1 .//2 .^(5:-1:0)
-numtraj = Int(5e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, NON2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u) .- u₀.*exp(1.0*(p[1]+0.5*p[2]^2))) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON2:", m)
-
+sim = test_convergence(dts,ensemble_prob,COM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=u₀*exp(1.0*(p[1]+0.5*p[2]^2))
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("COM:", sim.𝒪est[:weak_final])
 
 """
  Test Diagonal noise SDEs (iip)
@@ -250,64 +202,52 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.37
-
-println("RS1:", m)
-
-numtraj = Int(4e4)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("RS2:", m)
-
+sim = test_convergence(dts,ensemble_prob,RS1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.34 # order is 1.67
+println("RS1:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RS2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RS2:", sim.𝒪est[:weak_final])
 
 numtraj = Int(1e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, NON(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON:", m)
+sim = test_convergence(dts,ensemble_prob,NON(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,NON2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON2:", sim.𝒪est[:weak_final])
 
 numtraj = Int(6e4)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, COM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("COM:", m)
-
-
-numtraj = Int(1e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, NON2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(301//100)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON2:", m)
+sim = test_convergence(dts,ensemble_prob,COM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(301//100)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("COM:", sim.𝒪est[:weak_final])
diff --git a/test/weak_convergence/weak_strat_non_diagonal.jl b/test/weak_convergence/weak_strat_non_diagonal.jl
index 1512ad67c..952dfa58d 100644
--- a/test/weak_convergence/weak_strat_non_diagonal.jl
+++ b/test/weak_convergence/weak_strat_non_diagonal.jl
@@ -5,24 +5,14 @@
  NON, COM, NON2
 """
 
-
 import Statistics # for mean values of trajectories
 import LinearAlgebra # for the normn
 using StochasticDiffEq
 using Test
 using Random
+using DiffEqDevTools
 #using DiffEqGPU
 
-function generate_weak_solutions(prob, alg, dts, numtraj; ensemblealg=EnsembleThreads())
-  sols = []
-  for i in 1:length(dts)
-    sol = solve(prob,alg;ensemblealg=ensemblealg,dt=dts[i],adaptive=false,save_start=false,save_everystep=false,weak_timeseries_errors=false,weak_dense_errors=false,trajectories=Int(numtraj))
-    println(i)
-    push!(sols,sol)
-  end
-  return sols
-end
-
 function prob_func(prob, i, repeat)
     remake(prob,seed=seeds[i])
 end
@@ -60,64 +50,48 @@ seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, RS1(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("RS1:", m)
-
-
-numtraj = Int(1e6)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, RS2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("RS2:", m)
-
-numtraj = Int(1e6)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, NON(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
-m = log(errors[end]/errors[1])/log(dts[end]/dts[1])
-@test abs(m-2) < 0.3
-
-println("NON:", m)
-
-numtraj = Int(5e5)
-seed = 100
-Random.seed!(seed)
-seeds = rand(UInt, numtraj)
-
-_solutions = @time generate_weak_solutions(ensemble_prob, COM(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
-m = log(errors[end]/errors[2])/log(dts[end]/dts[2])
-@test abs(m-2) < 0.3 # tests are passing; problem might be not hard enough..
-
-println("COM:", m)
-
+sim = test_convergence(dts,ensemble_prob,RS1(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(2)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RS1:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,RS2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(2)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("RS2:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,NON(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(2)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON:", sim.𝒪est[:weak_final])
+
+sim = test_convergence(dts,ensemble_prob,NON2(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(2)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("NON2:", sim.𝒪est[:weak_final])
 
 numtraj = Int(5e5)
 seed = 100
 Random.seed!(seed)
 seeds = rand(UInt, numtraj)
 
-_solutions = @time generate_weak_solutions(ensemble_prob, NON2(), dts, numtraj, ensemblealg=EnsembleThreads())
-
-errors = [LinearAlgebra.norm(Statistics.mean(sol.u)-1//100*exp(2)) for sol in _solutions]
-m = log(errors[end]/errors[2])/log(dts[end]/dts[2])
-@test abs(m-2) < 0.3 # tests are passing; problem might be not hard enough..
-
-println("NON2:", m)
+sim = test_convergence(dts,ensemble_prob,COM(),
+    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
+    weak_timeseries_errors=false,weak_dense_errors=false,
+    expected_value=1//100*exp(2)
+    )
+@test abs(sim.𝒪est[:weak_final]-2) < 0.3
+println("COM:", sim.𝒪est[:weak_final])
+# COM tests are passing; problem might be not hard enough..