diff --git a/Project.toml b/Project.toml
index 82b0a555..c340d95f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
 Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
@@ -27,6 +28,7 @@ IntegralsCubaExt = "Cuba"
 IntegralsCubatureExt = "Cubature"
 IntegralsFastGaussQuadratureExt = "FastGaussQuadrature"
 IntegralsForwardDiffExt = "ForwardDiff"
+IntegralsGSLExt = "GSL"
 IntegralsZygoteExt = ["Zygote", "ChainRulesCore"]
 
 [compat]
@@ -40,6 +42,7 @@ Distributions = "0.25.87"
 FastGaussQuadrature = "0.5"
 FiniteDiff = "2.12"
 ForwardDiff = "0.10.19"
+GSL = "1"
 HCubature = "1.5"
 LinearAlgebra = "1.9"
 MonteCarloIntegration = "0.0.3, 0.1"
@@ -64,6 +67,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
@@ -72,4 +76,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Arblib", "SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"]
+test = ["Aqua", "Arblib", "SciMLSensitivity", "StaticArrays", "FiniteDiff", "GSL", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"]
diff --git a/ext/IntegralsGSLExt.jl b/ext/IntegralsGSLExt.jl
new file mode 100644
index 00000000..87babf6f
--- /dev/null
+++ b/ext/IntegralsGSLExt.jl
@@ -0,0 +1,48 @@
+module IntegralsGSLExt
+
+using GSL
+using Integrals
+using Integrals: IntegralCache
+
+mutable struct GSLCache{T}
+    value::T
+end
+getvalue(cache::GSLCache) = cache.value
+
+function Integrals.init_cacheval(alg::GSLIntegration{typeof(integration_cquad)}, prob::IntegralProblem)
+    ws = integration_cquad_workspace_alloc(alg.kws.wssize)
+    gslcache = GSLCache(ws)
+    finalizer(integration_cquad_workspace_free∘getvalue, gslcache)
+    result = Cdouble[0]
+    abserr = Cdouble[0]
+    nevals = C_NULL # Csize_t[0]
+    return (; gslcache, result, abserr, nevals)
+end
+
+function Integrals.__solvebp_call(cache::IntegralCache, alg::GSLIntegration{typeof(integration_cquad)}, sensealg, domain, p;
+    reltol = 1e-8, abstol = 1e-8, maxiters = nothing)
+
+    prob = Integrals.build_problem(cache)
+
+    if !all(isone∘length, domain)
+        error("GSLIntegration only accepts one-dimensional quadrature problems.")
+    end
+    @assert prob.f isa IntegralFunction
+
+    f = if isinplace(prob)
+        @assert isone(length(prob.f.integrand_prototype)) "GSL only supports scalar, real-valued integrands"
+        y = similar(prob.f.integrand_prototype, Cdouble)
+        x -> (prob.f(y, x, p); only(y))
+    else
+        x -> Cdouble(only(prob.f(x, p)))
+    end
+    # gslf = @gsl_function(f) # broken, see: https://github.com/JuliaMath/GSL.jl/pull/128
+    ptr = @cfunction($((x,p) -> f(x)), Cdouble, (Cdouble, Ptr{Cvoid}))
+    gslf = GC.@preserve ptr gsl_function(Base.unsafe_convert(Ptr{Cvoid},ptr), 0)
+    a, b = map(Cdouble∘only, domain)
+    (; gslcache, result, abserr, nevals) = cache.cacheval
+    integration_cquad(gslf, a, b, abstol, reltol, getvalue(gslcache), result, abserr, nevals)
+    return SciMLBase.build_solution(prob, alg, only(result), only(abserr), retcode = ReturnCode.Success)
+end
+
+end
diff --git a/src/Integrals.jl b/src/Integrals.jl
index 47faddc2..d0ee6c8b 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -221,5 +221,6 @@ export QuadGKJL, HCubatureJL, VEGAS, GaussLegendre, QuadratureRule, TrapezoidalR
 export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
 export CubatureJLh, CubatureJLp
 export ArblibJL
+export GSLIntegration
 
 end # module
diff --git a/src/algorithms.jl b/src/algorithms.jl
index d40dceaa..3c042c6a 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -370,3 +370,16 @@ end
 function ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL)
     return ArblibJL(check_analytic, take_prec, warn_on_no_convergence, opts)
 end
+
+
+"""
+    GSLIntegration(routine; kws...)
+
+One-dimensional quadrature of Float64-valued function using `routine` from GSL with
+additional arguments. For example `using Integrals, GSL; GSLIntegration(integration_cquad; wssize=100)`
+"""
+struct GSLIntegration{F,A<:NamedTuple} <: SciMLBase.AbstractIntegralAlgorithm
+    f::F
+    kws::A
+end
+GSLIntegration(f; kws...) = GSLIntegration(f, NamedTuple(kws))
diff --git a/test/interface_tests.jl b/test/interface_tests.jl
index 8bd9603a..ab0d5a30 100644
--- a/test/interface_tests.jl
+++ b/test/interface_tests.jl
@@ -1,5 +1,5 @@
 using Integrals
-using Cuba, Cubature, Arblib
+using Cuba, Cubature, Arblib, GSL
 using Test
 
 max_dim_test = 2
@@ -9,7 +9,7 @@ abstol = 1e-3
 
 
 algs = [QuadGKJL, HCubatureJL, CubatureJLh, CubatureJLp, #VEGAS, #CubaVegas,
-    CubaSUAVE, CubaDivonne, CubaCuhre]
+    CubaSUAVE, CubaDivonne, CubaCuhre, ArblibJL, () -> GSL(integration_cquad; wssize=100)]
 
 alg_req = Dict(QuadGKJL => (nout = 1, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true),
@@ -29,7 +29,8 @@ alg_req = Dict(QuadGKJL => (nout = 1, allows_batch = true, min_dim = 1, max_dim
         max_dim = Inf, allows_iip = true),
     CubaCuhre => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
         allows_iip = true),
-    ArblibJL => (nout=1, allows_batch=false, min_dim=1, max_dim=1, allows_iip=true))
+    ArblibJL => (nout=1, allows_batch=false, min_dim=1, max_dim=1, allows_iip=true,
+    GSLIntegration => (nout=1, allows_batch=false, min_dim=1, max_dim=1, allows_iip=true)))
 
 integrands = [
     (x, p) -> 1.0,