diff --git a/.gitignore b/.gitignore
index 1dfdc8a..9bd0ce0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,4 @@
 /Manifest.toml
-test.jl
-
 docs/Manifest.toml
+ext/Manifest.toml
 docs/build
-example_data/d_t2/
diff --git a/Project.toml b/Project.toml
index b025cac..8678e05 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,9 +17,11 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
 RipQP = "1e40b3f8-35eb-4cd8-8edd-3e515bb9de08"
+NonNegLeastSquares = "b7351bd1-99d9-5c5d-8786-f205a815c4d7"
 
 [extensions]
 gui_ext = "GLMakie"
+nnls_ext = "NonNegLeastSquares"
 tikhonov_jump_ext = "JuMP"
 tikhonov_ripqp_ext = ["RipQP", "QuadraticModels"]
 
@@ -28,6 +30,7 @@ DelimitedFiles = "1"
 GLMakie = "0.11, 0.12, 0.13"
 JuMP = "1"
 NativeFileDialog = "0.2"
+NonNegLeastSquares = "0.4.1"
 Optim = "1.10.0"
 PolygonOps = "0.1"
 QuadraticModels = "0.9"
@@ -43,11 +46,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = [
-    "Test",
-    "GLMakie",
-    "NMRInversions", 
-    "Random", 
-    "SparseArrays",
-    "LinearAlgebra"
-]
+test = ["Test", "GLMakie", "NMRInversions", "Random", "SparseArrays", "LinearAlgebra"]
diff --git a/docs/src/types_structs.md b/docs/src/types_structs.md
index 3cf08b2..df8fc0a 100644
--- a/docs/src/types_structs.md
+++ b/docs/src/types_structs.md
@@ -35,6 +35,14 @@ cdL1
 optim_nnls
 ```
 
+The following are also implemented through
+external package extensions.
+
+```@docs
+nnls
+jump_nnls
+```
+
 ## Finding optimal alpha
 These are methods for finding the optimal regularization parameter. 
 They can be used as input to the `invert` function as the 'alpha' argument
diff --git a/ext/Manifest.toml b/ext/Manifest.toml
index cde1d76..a449c2a 100644
--- a/ext/Manifest.toml
+++ b/ext/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.12.4"
+julia_version = "1.12.5"
 manifest_format = "2.0"
-project_hash = "36034e3e9c89d48d9c539f2b0060983a7836f13f"
+project_hash = "55fd022efc13f657a37564ea8957b17208f5633d"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a"
@@ -1069,16 +1069,18 @@ version = "7.10.0"
 deps = ["DelimitedFiles", "LinearAlgebra", "NativeFileDialog", "Optim", "Pkg", "PolygonOps", "SparseArrays"]
 path = ".."
 uuid = "55c20db2-0166-4687-95c3-62a9c7afb29b"
-version = "1.2.0"
+version = "1.3.0"
 
     [deps.NMRInversions.extensions]
     gui_ext = "GLMakie"
+    nnls_ext = "NonNegLeastSquares"
     tikhonov_jump_ext = "JuMP"
     tikhonov_ripqp_ext = ["RipQP", "QuadraticModels"]
 
     [deps.NMRInversions.weakdeps]
     GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
     JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+    NonNegLeastSquares = "b7351bd1-99d9-5c5d-8786-f205a815c4d7"
     QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
     RipQP = "1e40b3f8-35eb-4cd8-8edd-3e515bb9de08"
 
@@ -1110,6 +1112,12 @@ version = "1.1.1"
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.3.0"
 
+[[deps.NonNegLeastSquares]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "cdc11138e74a0dd0b82e7d64eb1350fdf049d3b1"
+uuid = "b7351bd1-99d9-5c5d-8786-f205a815c4d7"
+version = "0.4.1"
+
 [[deps.Observables]]
 git-tree-sha1 = "7438a59546cf62428fc9d1bc94729146d37a7225"
 uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
diff --git a/ext/Project.toml b/ext/Project.toml
index 02c05ec..7080efe 100644
--- a/ext/Project.toml
+++ b/ext/Project.toml
@@ -3,4 +3,5 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NMRInversions = "55c20db2-0166-4687-95c3-62a9c7afb29b"
 NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
+NonNegLeastSquares = "b7351bd1-99d9-5c5d-8786-f205a815c4d7"
 PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
diff --git a/ext/nnls_ext.jl b/ext/nnls_ext.jl
new file mode 100644
index 0000000..ce82ecb
--- /dev/null
+++ b/ext/nnls_ext.jl
@@ -0,0 +1,28 @@
+module nnls_ext
+
+using NMRInversions
+using LinearAlgebra, NonNegLeastSquares
+
+function NMRInversions.solve_regularization(K::AbstractMatrix, g::AbstractVector, α::Real, solver::nnls)
+
+    L = if isa(solver.L,Int)
+        NMRInversions.Γ(size(K, 2), solver.L)
+    else
+        if size(solver.L, 2) == size(K, 2) 
+            solver.L 
+        else
+            throw("Size mismatch between `solver.L` and the Kernel.")
+        end
+    end
+
+    A = [K; √(α) .* L]
+    b = [g; zeros(size(A, 1) - size(g, 1))]
+
+    f = vec(nonneg_lsq(A, b ; alg=solver.algorithm))
+
+    r = K * f - g
+
+    return f, r
+end
+
+end
diff --git a/ext/tikhonov_jump_ext.jl b/ext/tikhonov_jump_ext.jl
index 6bee666..d42c0d3 100644
--- a/ext/tikhonov_jump_ext.jl
+++ b/ext/tikhonov_jump_ext.jl
@@ -8,31 +8,22 @@ using NMRInversions, JuMP, LinearAlgebra, SparseArrays
 # import JuMP: @constraints
 # import JuMP: @objective
 
-export jump_nnls
-"""
-    jump_nnls(order, solver)
-Jump non-negative least squares method for tikhonov (L2) regularization, 
-implemented using the JuMP extension.
-All around effective, but can be slow for large problems, such as 2D inversions, 
-unless a powerful solver like gurobi is used.
-
-- `solver` is passed directly into JuMP. 
-- `order` determines the tikhonov finite-difference matrix. \
-If 0 is chosen, the identity matrix is used.
-
-This one is still experimental and not well-tested, 
-please submit an issue if you encounter any difficulties.
-"""
-struct jump_nnls <: regularization_solver 
-    order::Int
-    solver::Symbol
-end
 
 function NMRInversions.solve_regularization(K::AbstractMatrix, g::AbstractVector,
         α::Real, solver::Type{jump_nnls})
 
-    A = sparse([K; √(α) .* NMRInversions.Γ(size(K, 2), solver.order)])
-    b = sparse([g; zeros(size(A, 1) - size(g, 1))])
+    L = if isa(solver.L,Int)
+        NMRInversions.Γ(size(K, 2), solver.L)
+    else
+        if size(solver.L, 2) == size(K, 2) 
+            solver.L 
+        else
+            throw("Size mismatch between `solver.L` and the Kernel.")
+        end
+    end
+
+    A = [K; √(α) .* L]
+    b = [g; zeros(size(A, 1) - size(g, 1))]
 
     @eval model = JuMP.Model($(solver.solver).Optimizer)
     set_silent(model)
diff --git a/src/brd.jl b/src/brd.jl
index 903762d..1208eb7 100644
--- a/src/brd.jl
+++ b/src/brd.jl
@@ -1,7 +1,7 @@
 using Optim
 
 """
-    brd
+    brd()
 Solver for tikhonov (L2) regularization, following \
 [this paper](https://doi.org/10.1109/78.995059)
 [Venka2002](@cite).
diff --git a/src/optim_least_squares.jl b/src/optim_least_squares.jl
index 1c15499..af84379 100644
--- a/src/optim_least_squares.jl
+++ b/src/optim_least_squares.jl
@@ -2,7 +2,7 @@ using Optim
 
 export optim_nnls
 """
-    optim_nnls(order)
+    optim_nnls(; L, algorithm, opts)
 Simple non-negative least squares method for regularization, 
 implemented using Optim.jl.
 All around effective, but can be slow for large problems, such as 2D inversions.
@@ -31,7 +31,7 @@ function solve_regularization(K::AbstractMatrix, g::AbstractVector, α::Real, so
         if size(solver.L, 2) == size(K, 2) 
             solver.L 
         else
-            throw("Size mismatch between `optim_nnls.L` and the Kernel.")
+            throw("Size mismatch between `solver.L` and the Kernel.")
         end
     end
 
diff --git a/src/types.jl b/src/types.jl
index ba61556..b710f3f 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -191,3 +191,55 @@ mutable struct expfit_struct
     eqn::String
     title::String
 end
+
+
+# define structs for extension solvers
+function solve_regularization(K, g, α, solver::regularization_solver)
+    error("The package needed for $(typeof(solver)) is not loaded.")
+end
+
+export nnls
+"""
+    nnls(; L, algorithm)
+Non-negative least squares method for regularization, 
+implemented using the `NonNegLeastSquares` package extension.
+Decent speed for small, 1D regularizations, but very slow for larger, 2D problems.
+It can be used as a "solver" for invert function.
+
+- `L` determines the tikhonov matrix.\
+Can be either a matrix or an integer `n`. The latter will create a finite difference matrix\
+of order `n`, which means that the penalty term will be the `n`'th derivative of the distribution.\
+Defaults to `0`, where `L` becomes the Identity matrix.
+- `algorithm` is one of:
+    * `:nnls`
+    * `:fnnls`
+    * `:pivot`
+"""
+struct nnls <: regularization_solver 
+    L::Union{Int, AbstractMatrix}
+    algorithm::Symbol
+end
+nnls(;L=0, algorithm=:nnls) = nnls(L, algorithm)
+
+
+export jump_nnls
+"""
+    jump_nnls(L, solver)
+Jump non-negative least squares method for tikhonov (L2) regularization, 
+implemented using the JuMP package extension.
+All around effective, but can be slow for large problems, such as 2D inversions, 
+unless a powerful solver like gurobi is used.
+
+- `L` determines the tikhonov matrix.\
+Can be either a matrix or an integer `n`. The latter will create a finite difference matrix\
+of order `n`, which means that the penalty term will be the `n`'th derivative of the distribution.\
+Defaults to `0`, where `L` becomes the Identity matrix.
+- `solver` is passed directly into JuMP, see its documentation for available options. 
+
+This one is still experimental and not well-tested, 
+please submit an issue if you encounter any difficulties.
+"""
+struct jump_nnls <: regularization_solver 
+    order::Int
+    solver::Symbol
+end