diff --git a/Project.toml b/Project.toml
index 7cee55dea..6c129234a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 ConsoleProgressMonitor = "88cd18e8-d9cc-4ea6-8889-5259c0d15c8b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
@@ -24,6 +25,7 @@ ADTypes = "0.2.5"
 ArrayInterface = "7.6"
 ConsoleProgressMonitor = "0.1.1"
 DocStringExtensions = "0.9"
+LBFGSB = "0.4.1"
 LinearAlgebra = "1.10"
 Logging = "1.10"
 LoggingExtras = "0.4, 1"
diff --git a/src/Optimization.jl b/src/Optimization.jl
index a5e83efe7..7945a4c93 100644
--- a/src/Optimization.jl
+++ b/src/Optimization.jl
@@ -23,6 +23,7 @@ export ObjSense, MaxSense, MinSense
 
 include("utils.jl")
 include("state.jl")
+include("lbfgsb.jl")
 
 export solve
 
diff --git a/src/lbfgsb.jl b/src/lbfgsb.jl
new file mode 100644
index 000000000..8433e5a30
--- /dev/null
+++ b/src/lbfgsb.jl
@@ -0,0 +1,90 @@
+using Optimization.SciMLBase, LBFGSB
+
+"""
+$(TYPEDEF)
+
+[L-BFGS-B](https://en.wikipedia.org/wiki/Limited-memory_BFGS#L-BFGS-B) Nonlinear Optimization Code from [LBFGSB](https://github.com/Gnimuc/LBFGSB.jl/tree/master).
+It is a quasi-Newton optimization algorithm that supports bounds.
+
+References
+
+- R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208.
+- C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550 - 560.
+- J.L. Morales and J. Nocedal. L-BFGS-B: Remark on Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization (2011), to appear in ACM Transactions on Mathematical Software.
+"""
+@kwdef struct LBFGS
+    m::Int = 10
+end
+
+SciMLBase.supports_opt_cache_interface(::LBFGS) = true
+SciMLBase.allowsbounds(::LBFGS) = true
+# SciMLBase.requiresgradient(::LBFGS) = true
+
+function SciMLBase.__init(prob::SciMLBase.OptimizationProblem,
+        opt::LBFGS,
+        data = Optimization.DEFAULT_DATA; save_best = true,
+        callback = (args...) -> (false),
+        progress = false, kwargs...)
+    return OptimizationCache(prob, opt, data; save_best, callback, progress,
+        kwargs...)
+end
+
+function SciMLBase.__solve(cache::OptimizationCache{
+        F,
+        RC,
+        LB,
+        UB,
+        LC,
+        UC,
+        S,
+        O,
+        D,
+        P,
+        C
+}) where {
+        F,
+        RC,
+        LB,
+        UB,
+        LC,
+        UC,
+        S,
+        O <:
+        LBFGS,
+        D,
+        P,
+        C
+}
+    if cache.data != Optimization.DEFAULT_DATA
+        maxiters = length(cache.data)
+        data = cache.data
+    else
+        maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters)
+        data = Optimization.take(cache.data, maxiters)
+    end
+
+    local x
+
+    _loss = function (θ)
+        x = cache.f(θ, cache.p)
+        opt_state = Optimization.OptimizationState(u = θ, objective = x[1])
+        if cache.callback(opt_state, x...)
+            error("Optimization halted by callback.")
+        end
+        return x[1]
+    end
+
+    t0 = time()
+    if cache.lb !== nothing && cache.ub !== nothing
+        res = lbfgsb(_loss, cache.f.grad, cache.u0; m = cache.opt.m, maxiter = maxiters,
+            lb = cache.lb, ub = cache.ub)
+    else
+        res = lbfgsb(_loss, cache.f.grad, cache.u0; m = cache.opt.m, maxiter = maxiters)
+    end
+
+    t1 = time()
+    stats = Optimization.OptimizationStats(; iterations = maxiters,
+        time = t1 - t0, fevals = maxiters, gevals = maxiters)
+
+    return SciMLBase.build_solution(cache, cache.opt, res[2], res[1], stats = stats)
+end
diff --git a/test/lbfgsb.jl b/test/lbfgsb.jl
new file mode 100644
index 000000000..e1a471160
--- /dev/null
+++ b/test/lbfgsb.jl
@@ -0,0 +1,19 @@
+using Optimization
+using ForwardDiff, Zygote, ReverseDiff, FiniteDiff, Tracker
+using ModelingToolkit, Enzyme, Random
+
+x0 = zeros(2)
+rosenbrock(x, p = nothing) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
+l1 = rosenbrock(x0)
+
+optf = OptimizationFunction(rosenbrock, AutoForwardDiff())
+prob = OptimizationProblem(optf, x0)
+res = solve(prob, Optimization.LBFGS(), maxiters = 100)
+
+@test res.u≈[1.0, 1.0] atol=1e-3
+
+optf = OptimizationFunction(rosenbrock, AutoZygote())
+prob = OptimizationProblem(optf, x0, lb = [0.0, 0.0], ub = [0.3, 0.3])
+res = solve(prob, Optimization.LBFGS(), maxiters = 100)
+
+@test res.u≈[0.3, 0.09] atol=1e-3