diff --git a/Project.toml b/Project.toml
index acb3112ba..6fb50e1ec 100644
--- a/Project.toml
+++ b/Project.toml
@@ -73,7 +73,7 @@ Flux = "0.14.22, 0.15, 0.16"
 ForwardDiff = "0.10.36, 1"
 Functors = "0.4.12, 0.5"
 Hwloc = "3.3.0"
-Integrals = "4.5"
+Integrals = "4.5 - 4.6, 4.9.1"
 InteractiveUtils = "<0.0.1, 1"
 IntervalSets = "0.7.10"
 LineSearches = "7.3"
@@ -90,9 +90,9 @@ ModelingToolkit = "9.46, 10"
 MonteCarloMeasurements = "1.1"
 NeuralOperators = "0.5, 0.6"
 Optimisers = "0.3.3, 0.4"
-Optimization = "4, 5"
+Optimization = "5.1"
 OptimizationOptimJL = "0.4"
-OptimizationOptimisers = "0.3"
+OptimizationOptimisers = "0.3.14"
 OrdinaryDiffEq = "6.87"
 Printf = "1.10"
 QuasiMonteCarlo = "0.3.2"
@@ -110,7 +110,7 @@ Symbolics = "6.14"
 TensorBoardLogger = "0.1.24"
 Test = "1.10"
 WeightInitializers = "1.0.3"
-Zygote = "0.6.71, 0.7"
+Zygote = "0.6.71, 0.7.0 - 0.7.9"
 julia = "1.10"
 
 [extras]
diff --git a/test/NNODE_tests.jl b/test/NNODE_tests.jl
index a015df690..724c36ea9 100644
--- a/test/NNODE_tests.jl
+++ b/test/NNODE_tests.jl
@@ -294,7 +294,7 @@ end
 end
 
 @testitem "ODE Complex Numbers" tags = [:nnode] begin
-    using OrdinaryDiffEq, Random, Lux, Optimisers
+    using OrdinaryDiffEq, Random, Lux, Optimisers, Integrals
 
     Random.seed!(100)
 
@@ -338,7 +338,8 @@ end
         @test sol.u ≈ ground_truth.u rtol = 1.0e-1
     end
 
-    alg = NNODE(chain, Adam(0.01); strategy = QuadratureTraining())
+    # Use HCubatureJL explicitly to avoid Julia base bug with @cfunction
+    alg = NNODE(chain, Adam(0.01); strategy = QuadratureTraining(quadrature_alg = HCubatureJL()))
     @test_throws ErrorException solve(
         problem, alg; verbose = false, maxiters = 5000, saveat = 0.01
     )
diff --git a/test/NNPDE_tests.jl b/test/NNPDE_tests.jl
index adaa4636b..150209560 100644
--- a/test/NNPDE_tests.jl
+++ b/test/NNPDE_tests.jl
@@ -1,6 +1,6 @@
 @testsetup module NNPDE1TestSetup
 
-using NeuralPDE, Cubature, Integrals, QuasiMonteCarlo
+using NeuralPDE, Integrals, QuasiMonteCarlo
 
 # DataGen is Real: https://github.com/SciML/NeuralPDE.jl/issues/906
 @parameters x
@@ -18,8 +18,10 @@ function callback(p, l)
 end
 
 grid_strategy = GridTraining(0.1)
+# Use HCubatureJL (pure Julia) instead of CubatureJLh to avoid Julia base bug
+# with @cfunction type parameters (spvals typo in base/meta.jl partially_inline!)
 quadrature_strategy = QuadratureTraining(
-    quadrature_alg = CubatureJLh(),
+    quadrature_alg = HCubatureJL(),
     reltol = 1.0e3, abstol = 1.0e-3, maxiters = 50, batch = 100
 )
 stochastic_strategy = StochasticTraining(100; bcs_points = 50)
@@ -45,7 +47,7 @@ export callback, strategies
 end
 
 @testitem "Test Heterogeneous ODE" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Cubature, Integrals, QuasiMonteCarlo, DomainSets, Lux, Random, Optimisers
+    using Integrals, QuasiMonteCarlo, DomainSets, Lux, Random, Optimisers
 
     function simple_1d_ode(strategy)
         @parameters θ
@@ -118,8 +120,10 @@ end
     ]
 
     grid_strategy = GridTraining(0.1)
+    # Note: quadrature_strategy defined but not used in this test
+    # Using HCubatureJL instead of CubatureJLh to avoid Julia base bug
     quadrature_strategy = QuadratureTraining(
-        quadrature_alg = CubatureJLh(),
+        quadrature_alg = HCubatureJL(),
         reltol = 1.0e-3, abstol = 1.0e-3, maxiters = 50, batch = 100
     )
 
@@ -168,7 +172,7 @@ end
 end
 
 @testitem "PDE II: 2D Poisson" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
 
     function test_2d_poisson_equation(chain, strategy)
@@ -214,7 +218,9 @@ end
         test_2d_poisson_equation(chain, strategy)
     end
 
-    algs = [CubatureJLp()]
+    # Use HCubatureJL instead of CubatureJLp to avoid Julia base bug
+    # with @cfunction type parameters (spvals typo in base/meta.jl)
+    algs = [HCubatureJL()]
     @testset "$(nameof(typeof(alg)))" for alg in algs
         strategy = QuadratureTraining(
             quadrature_alg = alg, reltol = 1.0e-4,
@@ -225,7 +231,7 @@ end
 end
 
 @testitem "PDE III: 3rd-order ODE" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
     import OptimizationOptimJL: BFGS
 
@@ -297,7 +303,7 @@ end
 end
 
 @testitem "PDE IV: System of PDEs" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
 
     @parameters x, y
@@ -324,8 +330,10 @@ end
     chain1 = Chain(Dense(2, 15, tanh), Dense(15, 1))
     chain2 = Chain(Dense(2, 15, tanh), Dense(15, 1))
 
+    # Use HCubatureJL instead of CubatureJLh to avoid Julia base bug
+    # with @cfunction type parameters (spvals typo in base/meta.jl)
     quadrature_strategy = QuadratureTraining(
-        quadrature_alg = CubatureJLh(),
+        quadrature_alg = HCubatureJL(),
         reltol = 1.0e-3, abstol = 1.0e-3, maxiters = 50, batch = 100
     )
     chain = [chain1, chain2]
@@ -354,7 +362,7 @@ end
 end
 
 @testitem "PDE V: 2D Wave Equation" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals,
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals,
         LineSearches, Integrals
     import DomainSets: Interval, infimum, supremum
     import OptimizationOptimJL: BFGS
@@ -389,8 +397,10 @@ end
     phi = NeuralPDE.Phi(chain)
     derivative = NeuralPDE.numeric_derivative
 
+    # Use HCubatureJL instead of CubatureJLh to avoid Julia base bug
+    # with @cfunction type parameters (spvals typo in base/meta.jl)
     quadrature_strategy = QuadratureTraining(
-        quadrature_alg = CubatureJLh(),
+        quadrature_alg = HCubatureJL(),
         reltol = 1.0e-3, abstol = 1.0e-3, maxiters = 50, batch = 100
     )
 
@@ -423,7 +433,7 @@ end
 end
 
 @testitem "PDE VI: PDE with mixed derivative" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
     using OptimizationOptimJL: BFGS
     using LineSearches: BackTracking
@@ -467,7 +477,7 @@ end
 end
 
 @testitem "NNPDE: Translating from Flux" tags = [:nnpde1] setup = [NNPDE1TestSetup] begin
-    using Lux, Random, Optimisers, DomainSets, Cubature, QuasiMonteCarlo, Integrals
+    using Lux, Random, Optimisers, DomainSets, QuasiMonteCarlo, Integrals
     import DomainSets: Interval, infimum, supremum
     import OptimizationOptimJL: BFGS
     import Flux
@@ -481,7 +491,8 @@ end
     domains = [θ ∈ Interval(0.0, 1.0)]
 
     chain = Flux.Chain(Flux.Dense(1, 12, Flux.σ), Flux.Dense(12, 1))
-    discretization = PhysicsInformedNN(chain, QuadratureTraining())
+    # Use HCubatureJL explicitly to avoid Julia base bug with @cfunction
+    discretization = PhysicsInformedNN(chain, QuadratureTraining(quadrature_alg = HCubatureJL()))
     @test discretization.chain isa Lux.AbstractLuxLayer
 
     @named pde_system = PDESystem(eq, bcs, domains, [θ], [u])
diff --git a/test/direct_function_tests.jl b/test/direct_function_tests.jl
index fef14d5e6..78f537c0d 100644
--- a/test/direct_function_tests.jl
+++ b/test/direct_function_tests.jl
@@ -138,7 +138,7 @@ end
 end
 
 @testitem "Empty boundary condition [] fails in solve phase" begin
-    using NeuralPDE, Optimization, OptimizationOptimisers, Lux
+    using NeuralPDE, Optimization, OptimizationOptimisers, Lux, Integrals
     @parameters x
     @variables u(..)
 
@@ -151,7 +151,8 @@ end
             GridTraining(0.01),
             StochasticTraining(1000),
             QuasiRandomTraining(1000),
-            QuadratureTraining(),
+            # Use HCubatureJL explicitly to avoid Julia base bug with @cfunction
+            QuadratureTraining(quadrature_alg = HCubatureJL()),
         )
         discretization = PhysicsInformedNN(chain, strategy)
         @named pde_system = PDESystem(eq, bc, domain, [x], [u(x)])
diff --git a/test/neural_adapter_tests.jl b/test/neural_adapter_tests.jl
index 219fd1d59..a7c86ae0c 100644
--- a/test/neural_adapter_tests.jl
+++ b/test/neural_adapter_tests.jl
@@ -12,7 +12,7 @@ end
 
 @testitem "Neural Adapter: 2D Poisson" tags = [:neuraladapter] setup = [NeuralAdapterTestSetup] begin
     using Optimization, Lux, OptimizationOptimisers, Statistics, ComponentArrays, Random,
-        LinearAlgebra
+        LinearAlgebra, Integrals
     import DomainSets: Interval, infimum, supremum
 
     Random.seed!(100)
@@ -34,7 +34,9 @@ end
     ]
     # Space and time domains
     domains = [x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]
+    # Use HCubatureJL explicitly to avoid Julia base bug with @cfunction
     quadrature_strategy = QuadratureTraining(
+        quadrature_alg = HCubatureJL(),
         reltol = 1.0e-3, abstol = 1.0e-6, maxiters = 50, batch = 100
     )
     chain1 = Chain(Dense(2, 8, tanh), Dense(8, 8, tanh), Dense(8, 1))
@@ -60,7 +62,9 @@ end
     loss(cord, θ) = first(chain2(cord, θ, st)) .- phi(cord, res.u)
 
     grid_strategy = GridTraining(0.05)
+    # Use HCubatureJL explicitly to avoid Julia base bug with @cfunction
     quadrature_strategy = QuadratureTraining(
+        quadrature_alg = HCubatureJL(),
         reltol = 1.0e-3, abstol = 1.0e-6, maxiters = 50, batch = 100
     )
     stochastic_strategy = StochasticTraining(1000)