diff --git a/README.md b/README.md
index e8155a7..f934eac 100644
--- a/README.md
+++ b/README.md
@@ -34,7 +34,8 @@ into the problem as the fourth argument of `IntegralProblem`.
 using Integrals
 f(x, p) = sin(x * p)
 p = 1.7
-prob = IntegralProblem(f, -2, 5, p)
+domain = (-2, 5) # (lb, ub)
+prob = IntegralProblem(f, domain, p)
 sol = solve(prob, QuadGKJL())
 ```
 
@@ -45,7 +46,8 @@ the bounds giving the interval of integration for `x[i]`.
 ```julia
 using Integrals
 f(x, p) = sum(sin.(x))
-prob = IntegralProblem(f, ones(2), 3ones(2))
+domain = (ones(2), 3ones(2)) # (lb, ub)
+prob = IntegralProblem(f, domain)
 sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
 ```
 
@@ -60,7 +62,8 @@ function f(dx, x, p)
         dx[i] = sum(sin, @view(x[:, i]))
     end
 end
-prob = IntegralProblem(BatchIntegralFunction(f, zeros(0)), ones(2), 3ones(2))
+domain = (ones(2), 3ones(2)) # (lb, ub)
+prob = IntegralProblem(BatchIntegralFunction(f, zeros(0)), domain)
 sol = solve(prob, CubatureJLh(), reltol = 1e-3, abstol = 1e-3)
 ```
 
@@ -68,6 +71,6 @@ If we would like to compare the results against Cuba.jl's `Cuhre` method, then
 the change is a one-argument change:
 
 ```julia
-using IntegralsCuba
+using Cuba
 sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3)
 ```
diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md
index e7fc0a1..fd2aa84 100644
--- a/docs/src/basics/SampledIntegralProblem.md
+++ b/docs/src/basics/SampledIntegralProblem.md
@@ -28,6 +28,11 @@ The exact answer is of course \$ 1/3 \$.
 
 ## Details
 
+```@docs
+SciMLBase.SampledIntegralProblem
+solve(::SampledIntegralProblem, ::SciMLBase.AbstractIntegralAlgorithm)
+```
+
 ### Non-equidistant grids
 
 If the sampling points `x` are provided as an `AbstractRange`
diff --git a/docs/src/solvers/IntegralSolvers.md b/docs/src/solvers/IntegralSolvers.md
index 85f5017..17df07f 100644
--- a/docs/src/solvers/IntegralSolvers.md
+++ b/docs/src/solvers/IntegralSolvers.md
@@ -21,6 +21,8 @@ The following algorithms are available:
 ```@docs
 QuadGKJL
 HCubatureJL
+CubatureJLp
+CubatureJLh
 VEGAS
 VEGASMC
 CubaVegas
diff --git a/docs/src/tutorials/caching_interface.md b/docs/src/tutorials/caching_interface.md
index 61e5f5c..415e536 100644
--- a/docs/src/tutorials/caching_interface.md
+++ b/docs/src/tutorials/caching_interface.md
@@ -6,7 +6,7 @@ problems. For example, if one is going to solve the same integral for several pa
 ```julia
 using Integrals
 
-prob = IntegralProblem((x, p) -> sin(x * p), 0, 1, 14.0)
+prob = IntegralProblem((x, p) -> sin(x * p), (0, 1), 14.0)
 alg = QuadGKJL()
 
 solve(prob, alg)
@@ -33,7 +33,7 @@ This uses the [SciML `init` interface](https://docs.sciml.ai/SciMLBase/stable/in
 ```@example cache1
 using Integrals
 
-prob = IntegralProblem((x, p) -> sin(x * p), 0, 1, 14.0)
+prob = IntegralProblem((x, p) -> sin(x * p), (0, 1), 14.0)
 alg = QuadGKJL()
 
 cache = init(prob, alg)
@@ -45,7 +45,7 @@ cache.p = 15.0
 sol2 = solve!(cache)
 ```
 
-The caching interface is intended for updating `p`, `lb`, `ub`, `nout`, and `batch`.
+The caching interface is intended for updating `p` and `domain`.
 Note that the types of these variables is not allowed to change.
 If it is necessary to change the integrand `f` instead of defining a new
 `IntegralProblem`, consider using
diff --git a/docs/src/tutorials/differentiating_integrals.md b/docs/src/tutorials/differentiating_integrals.md
index 413d787..763f8a0 100644
--- a/docs/src/tutorials/differentiating_integrals.md
+++ b/docs/src/tutorials/differentiating_integrals.md
@@ -8,17 +8,16 @@ and Zygote.jl for reverse-mode automatic differentiation. For example:
 ```@example AD
 using Integrals, ForwardDiff, FiniteDiff, Zygote, Cuba
 f(x, p) = sum(sin.(x .* p))
-lb = ones(2)
-ub = 3ones(2)
+domain = (ones(2), 3ones(2)) # (lb, ub)
 p = ones(2)
 
 function testf(p)
-    prob = IntegralProblem(f, lb, ub, p)
+    prob = IntegralProblem(f, domain, p)
     sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1])
 end
 testf(p)
-#dp1 = Zygote.gradient(testf,p) #broken
+dp1 = Zygote.gradient(testf,p)
 dp2 = FiniteDiff.finite_difference_gradient(testf, p)
 dp3 = ForwardDiff.gradient(testf, p)
-dp2 ≈ dp3 #dp1[1] ≈ dp2 ≈ dp3
+dp1[1] ≈ dp2 ≈ dp3
 ```
diff --git a/docs/src/tutorials/numerical_integrals.md b/docs/src/tutorials/numerical_integrals.md
index fe9cf90..35ccf42 100644
--- a/docs/src/tutorials/numerical_integrals.md
+++ b/docs/src/tutorials/numerical_integrals.md
@@ -12,7 +12,8 @@ We can numerically approximate this integral:
 ```@example integrate1
 using Integrals
 f(u, p) = sum(sin.(u))
-prob = IntegralProblem(f, ones(3), 3ones(3))
+domain = (ones(3), 3ones(3)) # (lb, ub)
+prob = IntegralProblem(f, domain)
 sol = solve(prob, HCubatureJL(); reltol = 1e-3, abstol = 1e-3)
 sol.u
 ```
@@ -39,7 +40,8 @@ For example, we also want to evaluate:
 ```@example integrate2
 using Integrals
 f(u, p) = [sum(sin.(u)), sum(cos.(u))]
-prob = IntegralProblem(f, ones(3), 3ones(3))
+domain = (ones(3), 3ones(3)) # (lb, ub)
+prob = IntegralProblem(f, domain)
 sol = solve(prob, HCubatureJL(); reltol = 1e-3, abstol = 1e-3)
 sol.u
 ```
@@ -59,7 +61,8 @@ function f(y, u, p)
     y[2] = sum(cos.(u))
 end
 prototype = zeros(2)
-prob = IntegralProblem(IntegralFunction(f, prototype), ones(3), 3ones(3))
+domain = (ones(3), 3ones(3)) # (lb, ub)
+prob = IntegralProblem(IntegralFunction(f, prototype), domain)
 sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)
 sol.u
 ```
@@ -84,7 +87,8 @@ function f(y, u, p)
     end
 end
 prototype = zeros(2, 0)
-prob = IntegralProblem(BatchIntegralFunction(f, prototype), ones(3), 3ones(3))
+domain = (ones(3), 3ones(3)) # (lb, ub)
+prob = IntegralProblem(BatchIntegralFunction(f, prototype), domain)
 sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)
 sol.u
 ```
@@ -103,7 +107,8 @@ the change is a one-argument change:
 using Integrals
 using Cuba
 f(u, p) = sum(sin.(u))
-prob = IntegralProblem(f, ones(3), 3ones(3))
+domain = (ones(3), 3ones(3)) # (lb, ub)
+prob = IntegralProblem(f, domain)
 sol = solve(prob, CubaCuhre(); reltol = 1e-3, abstol = 1e-3)
 sol.u
 ```
@@ -118,7 +123,7 @@ For example, we can create our own sine function by integrating the cosine funct
 
 ```@example integrate6
 using Integrals
-my_sin(x) = solve(IntegralProblem((x, p) -> cos(x), 0.0, x), QuadGKJL()).u
+my_sin(x) = solve(IntegralProblem((x, p) -> cos(x), (0.0, x)), QuadGKJL()).u
 x = 0:0.1:(2 * pi)
 @. my_sin(x) ≈ sin(x)
 ```
@@ -152,7 +157,7 @@ using Distributions
 using Integrals
 dist = MvNormal(ones(2))
 f = (x, p) -> pdf(dist, x)
-(lb, ub) = ([-Inf, 0.0], [Inf, Inf])
-prob = IntegralProblem(f, lb, ub)
+domain = ([-Inf, 0.0], [Inf, Inf]) # (lb, ub)
+prob = IntegralProblem(f, domain)
 solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
 ```
diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl
index 111431e..e391c65 100644
--- a/ext/IntegralsCubaExt.jl
+++ b/ext/IntegralsCubaExt.jl
@@ -8,8 +8,8 @@ import Integrals: transformation_if_inf,
 function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
         sensealg,
         domain, p;
-        reltol = 1e-8, abstol = 1e-8,
-        maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int))
+        reltol = 1e-4, abstol = 1e-12,
+        maxiters = 1000000)
     @assert maxiters>=1000 "maxiters for $alg should be larger than 1000"
     lb, ub = domain
     mid = (lb + ub) / 2
@@ -19,16 +19,16 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
     # we could support other types by multiplying by the jacobian determinant at the end
 
     if prob.f isa BatchIntegralFunction
+        nvec = min(maxiters, prob.f.max_batch)
         # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when
         # batching
-        (nvec = prob.f.max_batch) > 1 ||
-            throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
+        nvec > 1 || throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
 
         if mid isa Real
-            _x = zeros(typeof(mid), prob.f.max_batch)
+            _x = zeros(typeof(mid), nvec)
             scale = x -> scale_x!(resize!(_x, length(x)), ub, lb, vec(x))
         else
-            _x = zeros(eltype(mid), length(mid), prob.f.max_batch)
+            _x = zeros(eltype(mid), length(mid), nvec)
             scale = x -> scale_x!(view(_x, :, 1:size(x, 2)), ub, lb, x)
         end
 
diff --git a/src/common.jl b/src/common.jl
index c1d4a03..cd31327 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -55,10 +55,12 @@ function Base.getproperty(cache::IntegralCache, name::Symbol)
 end
 function Base.setproperty!(cache::IntegralCache, name::Symbol, x)
     if name === :lb
+        @warn "updating lb is deprecated by replacing domain"
         lb, ub = cache.domain
         setfield!(cache, :domain, (oftype(lb, x), ub))
         return x
     elseif name === :ub
+        @warn "updating ub is deprecated by replacing domain"
         lb, ub = cache.domain
         setfield!(cache, :domain, (lb, oftype(ub, x)))
         return x
diff --git a/test/interface_tests.jl b/test/interface_tests.jl
index 1a0137c..19bd1c6 100644
--- a/test/interface_tests.jl
+++ b/test/interface_tests.jl
@@ -24,8 +24,8 @@ alg_req = Dict(
         max_dim = Inf, allows_iip = true),
     CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1,
         max_dim = Inf, allows_iip = true),
-    # CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
-    #     allows_iip = true),
+    CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
+        allows_iip = true),
     CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
         allows_iip = true),
     CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2,