diff --git a/Project.toml b/Project.toml
index 528576d..c41ed44 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,14 +23,14 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 Aqua = "0.8"
 BenchmarkTools = "1"
 ComponentArrays = "0.15"
+Cuba = "2"
+Cubature = "1.5"
 DiffEqBase = "6.100"
 DiffEqNoiseProcess = "5"
 Distributions = "0.23, 0.24, 0.25"
 FiniteDiff = "2"
 ForwardDiff = "0.10"
-Integrals = "3.1"
-IntegralsCuba = "0.2"
-IntegralsCubature = "0.2"
+Integrals = "4"
 KernelDensity = "0.6"
 LinearAlgebra = "1"
 OrdinaryDiffEq = "6"
@@ -56,11 +56,11 @@ julia = "1.6"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
-IntegralsCuba = "e00cd5f1-6337-4131-8b37-28b2fe4cd6cb"
-IntegralsCubature = "c31f79ba-6e32-46d4-a52f-182a8ac42a54"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
@@ -74,4 +74,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
 [targets]
-test = ["Test", "Aqua", "SafeTestsets", "TestExtras", "Pkg", "PkgBenchmark", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "OrdinaryDiffEq", "ComponentArrays", "ForwardDiff", "Integrals", "IntegralsCubature", "IntegralsCuba", "FiniteDiff", "RecursiveArrayTools", "StochasticDiffEq"]
+test = ["Test", "Aqua", "SafeTestsets", "TestExtras", "Pkg", "PkgBenchmark", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "OrdinaryDiffEq", "ComponentArrays", "ForwardDiff", "Integrals", "Cubature", "Cuba", "FiniteDiff", "RecursiveArrayTools", "StochasticDiffEq"]
diff --git a/docs/src/tutorials/gpu_bayesian.md b/docs/src/tutorials/gpu_bayesian.md
index 3b5af9f..39c3c11 100644
--- a/docs/src/tutorials/gpu_bayesian.md
+++ b/docs/src/tutorials/gpu_bayesian.md
@@ -22,7 +22,7 @@ Let's start by importing all the necessary libraries:
 using OrdinaryDiffEq
 using Turing, MCMCChains, Distributions
 using KernelDensity
-using SciMLExpectations, IntegralsCuba
+using SciMLExpectations, Cuba
 using DiffEqGPU
 using Plots, StatsPlots
 using Random;
diff --git a/docs/src/tutorials/introduction.md b/docs/src/tutorials/introduction.md
index 80be7ad..27f9291 100644
--- a/docs/src/tutorials/introduction.md
+++ b/docs/src/tutorials/introduction.md
@@ -367,7 +367,7 @@ It is also possible to solve the various simulations in parallel by using the `b
 The default quadrature algorithm to solve `ExpectationProblem` does not support batch-mode evaluation. So, we first load dependencies for additional quadrature algorithms
 
 ```@example introduction
-using IntegralsCuba
+using Cuba
 ```
 
 We then solve our expectation as before, using a `batch=10` multi-thread parallelization via `EnsembleThreads()` of Cuba's SUAVE algorithm. We also introduce additional uncertainty in the model parameter.
diff --git a/docs/src/tutorials/process_noise.md b/docs/src/tutorials/process_noise.md
index 714a215..b7c3696 100644
--- a/docs/src/tutorials/process_noise.md
+++ b/docs/src/tutorials/process_noise.md
@@ -5,7 +5,7 @@ This done by representing the Wiener process using the [Kosambi–Karhunen–Lo
 
 ```@example process_noise
 using SciMLExpectations
-using IntegralsCuba
+using Cuba
 using StochasticDiffEq
 using DiffEqNoiseProcess
 using Distributions
diff --git a/src/expectation.jl b/src/expectation.jl
index b21d5ba..edb6fde 100644
--- a/src/expectation.jl
+++ b/src/expectation.jl
@@ -226,7 +226,8 @@ function integrate(quadalg, adalg::AbstractExpectationADAlgorithm, f, lb::TB, ub
                    nout = 1, batch = 0,
                    kwargs...) where {TB}
     #TODO check batch iip type stability w/ IntegralProblem{XXXX}
-    prob = IntegralProblem{batch > 0}(f, lb, ub, p; nout = nout, batch = batch)
+    batch = batch==0 ? nothing : batch
+    prob = IntegralProblem(f, lb, ub, p; nout = nout, batch = batch)
     solve(prob, quadalg; kwargs...)
 end
 
diff --git a/src/system_utils.jl b/src/system_utils.jl
index 993a536..f1c60aa 100644
--- a/src/system_utils.jl
+++ b/src/system_utils.jl
@@ -26,12 +26,6 @@ function SystemMap(prob, alg::SciMLBase.AbstractODEAlgorithm,
     SystemMap(prob, alg, ensemblealg, kwargs)
 end
 
-function (sm::SystemMap{DT})(u0, p) where {DT}
-    prob::DT = remake(sm.prob,
-                      u0 = convert(typeof(sm.prob.u0), u0),
-                      p = convert(typeof(sm.prob.p), p))
-    solve(prob, sm.alg; sm.kwargs...)
-end
 
 """
 ```julia
diff --git a/test/processnoise.jl b/test/processnoise.jl
index a3375b2..5ad57f2 100644
--- a/test/processnoise.jl
+++ b/test/processnoise.jl
@@ -1,7 +1,7 @@
 using Test, TestExtras
 
 using SciMLExpectations
-using IntegralsCuba
+using Cuba
 using StochasticDiffEq
 using DiffEqNoiseProcess
 using Distributions
diff --git a/test/solve.jl b/test/solve.jl
index 10bc5f6..a80b4b0 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -1,6 +1,6 @@
 using Test, TestExtras,
       SciMLExpectations, OrdinaryDiffEq, Distributions,
-      Integrals, IntegralsCubature, IntegralsCuba,
+      Integrals, Cubature, Cuba,
       StaticArrays, ComponentArrays, Random
 
 quadalgs = [