diff --git a/Project.toml b/Project.toml
index d111c08..3425252 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,9 +4,11 @@ authors = ["QC-Devs"]
 version = "0.1.0"
 
 [deps]
+Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index e94421d..1242b61 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,18 +2,28 @@
 
 julia_version = "1.10.3"
 manifest_format = "2.0"
-project_hash = "c8d2fd08fdad193f86f7d8e90ce12b1cd2ada39a"
+project_hash = "3de629beb68b66b2fe028bd786bd890cb6ccca4a"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "016833eb52ba2d6bea9fcb50ca295980e728ee24"
+git-tree-sha1 = "daf26bbdec60d9ca1c0003b70f389d821ddb4224"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "0.2.7"
+version = "1.2.1"
+weakdeps = ["ChainRulesCore", "EnzymeCore"]
+
+    [deps.ADTypes.extensions]
+    ADTypesChainRulesCoreExt = "ChainRulesCore"
+    ADTypesEnzymeCoreExt = "EnzymeCore"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
 uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
 version = "0.0.1"
 
+[[deps.AbstractTrees]]
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
+uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+version = "0.4.5"
+
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown", "Test"]
 git-tree-sha1 = "c0d491ef0b135fd7d63cbc6404286bc633329425"
@@ -111,12 +121,28 @@ git-tree-sha1 = "585a387a490f1c4bd88be67eea15b93da5e85db7"
 uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
 version = "0.2.5"
 
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra"]
+git-tree-sha1 = "575cd02e080939a33b6df6c5853d14924c08e35b"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.23.0"
+weakdeps = ["SparseArrays"]
+
+    [deps.ChainRulesCore.extensions]
+    ChainRulesCoreSparseArraysExt = "SparseArrays"
+
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
 git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1"
 uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
 version = "0.1.12"
 
+[[deps.CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.4"
+
 [[deps.CommonSolve]]
 git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
 uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
@@ -199,9 +225,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"]
-git-tree-sha1 = "03b9555f4c3a7c2f530bb1ae13e85719c632f74e"
+git-tree-sha1 = "37d49a1f8eedfe68b7622075ff3ebe3dd0e8f327"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.151.1"
+version = "6.151.2"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -241,21 +267,57 @@ git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
 version = "1.15.1"
 
+[[deps.DifferentiationInterface]]
+deps = ["ADTypes", "Compat", "DocStringExtensions", "FillArrays", "LinearAlgebra", "PackageExtensionCompat", "SparseArrays", "SparseMatrixColorings"]
+git-tree-sha1 = "b4ce591b91f500126c9c5030aa5b8a710c7db7a0"
+uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+version = "0.4.2"
+
+    [deps.DifferentiationInterface.extensions]
+    DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
+    DifferentiationInterfaceDiffractorExt = "Diffractor"
+    DifferentiationInterfaceEnzymeExt = "Enzyme"
+    DifferentiationInterfaceFastDifferentiationExt = "FastDifferentiation"
+    DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
+    DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
+    DifferentiationInterfaceForwardDiffExt = "ForwardDiff"
+    DifferentiationInterfacePolyesterForwardDiffExt = "PolyesterForwardDiff"
+    DifferentiationInterfaceReverseDiffExt = "ReverseDiff"
+    DifferentiationInterfaceSymbolicsExt = "Symbolics"
+    DifferentiationInterfaceTapirExt = "Tapir"
+    DifferentiationInterfaceTrackerExt = "Tracker"
+    DifferentiationInterfaceZygoteExt = "Zygote"
+
+    [deps.DifferentiationInterface.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
+    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+    FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
+    FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+    FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+    Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+    Tapir = "07d77754-e150-4737-8c94-cd238a1fb45b"
+    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
-git-tree-sha1 = "c36550cb29cbe373e95b3f40486b9a4148f89ffd"
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.9.2"
+version = "0.9.3"
 
 [[deps.Documenter]]
-deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "6030186b00a38e9d0434518627426570aac2ef95"
+deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
+git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.23"
+version = "1.4.1"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -268,14 +330,26 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "1bc328eec34ffd80357f84a84bb30e4374e9bd60"
+git-tree-sha1 = "0910982db2490a20f81dc7db5d4bbea236c027b3"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.6.6"
+version = "0.7.3"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
     AdaptExt = "Adapt"
 
+[[deps.Evolutionary]]
+deps = ["LinearAlgebra", "NLSolversBase", "Random", "StackViews", "Statistics", "UnPack"]
+git-tree-sha1 = "495cc997e315ee71d2361bf4c90ef21dc12e8fd1"
+uuid = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
+version = "0.11.1"
+
+[[deps.Expat_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7"
+uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
+version = "2.6.2+0"
+
 [[deps.ExprTools]]
 git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
@@ -294,9 +368,9 @@ version = "0.3.2"
 
 [[deps.FastLapackInterface]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "f4102aab9c7df8691ed09f9c42e34f5ab5458ab9"
+git-tree-sha1 = "cbf5edddb61a43669710cbc2241bc08b36d9e660"
 uuid = "29a986be-02c6-4525-aec4-84b980013641"
-version = "2.0.3"
+version = "2.0.4"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
@@ -364,6 +438,18 @@ git-tree-sha1 = "ec632f177c0d990e64d955ccc1b8c04c485a0950"
 uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
 version = "0.1.6"
 
+[[deps.Git]]
+deps = ["Git_jll"]
+git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e"
+uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
+version = "1.3.1"
+
+[[deps.Git_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
+git-tree-sha1 = "d18fb8a1f3609361ebda9bf029b60fd0f120c809"
+uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
+version = "2.44.0+2"
+
 [[deps.Graphs]]
 deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
 git-tree-sha1 = "4f2b57488ac7ee16124396de4f2bbdd51b2602ad"
@@ -378,9 +464,9 @@ version = "0.1.16"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
-git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
+git-tree-sha1 = "8b72179abc660bfab5e28472e019392b97d0985c"
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
-version = "0.2.2"
+version = "0.2.4"
 
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
@@ -388,9 +474,9 @@ uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 version = "0.1.1"
 
 [[deps.Inflate]]
-git-tree-sha1 = "ea8031dea4aff6bd41f1df8f2fdfb25b33626381"
+git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d"
 uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
-version = "0.1.4"
+version = "0.1.5"
 
 [[deps.IntelOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -430,15 +516,15 @@ version = "1.5.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
+git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.3"
+version = "0.21.4"
 
 [[deps.KLU]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"]
-git-tree-sha1 = "884c2968c2e8e7e6bf5956af88cb46aa745c854b"
+git-tree-sha1 = "07649c499349dad9f08dde4243a4c597064663e9"
 uuid = "ef3ab10e-7fda-4108-b977-705223b18434"
-version = "0.4.1"
+version = "0.6.0"
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]
@@ -452,16 +538,29 @@ git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277"
 uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
 version = "0.1.15"
 
+[[deps.LazilyInitializedFields]]
+git-tree-sha1 = "8f7f3cabab0fd1800699663533b6d5cb3fc0e612"
+uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf"
+version = "1.2.2"
+
 [[deps.LazyArrays]]
-deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays"]
-git-tree-sha1 = "35079a6a869eecace778bcda8641f9a54ca3a828"
+deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "SparseArrays"]
+git-tree-sha1 = "1567f3b9c49a8249c0921a6c29c3caddecf77383"
 uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-version = "1.10.0"
-weakdeps = ["StaticArrays"]
+version = "2.0.2"
 
     [deps.LazyArrays.extensions]
+    LazyArraysBandedMatricesExt = "BandedMatrices"
+    LazyArraysBlockArraysExt = "BlockArrays"
+    LazyArraysBlockBandedMatricesExt = "BlockBandedMatrices"
     LazyArraysStaticArraysExt = "StaticArrays"
 
+    [deps.LazyArrays.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [[deps.LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
 uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
@@ -493,6 +592,12 @@ version = "1.11.0+1"
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[[deps.Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.17.0+0"
+
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
 git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d"
@@ -504,15 +609,16 @@ deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LinearSolve]]
-deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "StaticArraysCore", "UnPack"]
-git-tree-sha1 = "6f8e084deabe3189416c4e505b1c53e1b590cae8"
+deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "StaticArraysCore", "UnPack"]
+git-tree-sha1 = "7648cc20100504f4b453917aacc8520e9c0ecfb3"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "2.22.1"
+version = "2.30.1"
 
     [deps.LinearSolve.extensions]
     LinearSolveBandedMatricesExt = "BandedMatrices"
     LinearSolveBlockDiagonalsExt = "BlockDiagonals"
     LinearSolveCUDAExt = "CUDA"
+    LinearSolveCUDSSExt = "CUDSS"
     LinearSolveEnzymeExt = ["Enzyme", "EnzymeCore"]
     LinearSolveFastAlmostBandedMatricesExt = ["FastAlmostBandedMatrices"]
     LinearSolveHYPREExt = "HYPRE"
@@ -527,6 +633,7 @@ version = "2.22.1"
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
     EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
     FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
@@ -562,16 +669,12 @@ deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensio
 git-tree-sha1 = "8f6786d8b2b3248d79db3ad359ce95382d5a6df8"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 version = "0.12.170"
+weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"]
 
     [deps.LoopVectorization.extensions]
     ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]
     SpecialFunctionsExt = "SpecialFunctions"
 
-    [deps.LoopVectorization.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-    SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-
 [[deps.MKL_jll]]
 deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"]
 git-tree-sha1 = "80b2833b56d466b3858d565adcd16a4a05f2089b"
@@ -593,17 +696,17 @@ version = "0.1.8"
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
-[[deps.MatrixFactorizations]]
-deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"]
-git-tree-sha1 = "6731e0574fa5ee21c02733e397beb133df90de35"
-uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
-version = "2.2.0"
+[[deps.MarkdownAST]]
+deps = ["AbstractTrees", "Markdown"]
+git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899"
+uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391"
+version = "0.1.2"
 
 [[deps.MaybeInplace]]
 deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "SparseArrays"]
-git-tree-sha1 = "b1f2f92feb0bc201e91c155ef575bcc7d9cc3526"
+git-tree-sha1 = "1b9e613f2ca3b6cdcbfe36381e17ca2b66d4b3a1"
 uuid = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
-version = "0.1.2"
+version = "0.1.3"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -623,7 +726,7 @@ uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 version = "0.2.4"
 
 [[deps.NICE]]
-deps = ["ForwardDiff", "LinearAlgebra", "NonlinearSolve", "SciMLBase", "Test"]
+deps = ["Evolutionary", "ForwardDiff", "LinearAlgebra", "NonlinearSolve", "Reexport", "SciMLBase", "Test"]
 path = ".."
 uuid = "170674ed-ea2e-4c97-a263-9f057371600a"
 version = "0.1.0"
@@ -646,9 +749,9 @@ version = "1.2.0"
 
 [[deps.NonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "LazyArrays", "LineSearches", "LinearAlgebra", "LinearSolve", "MaybeInplace", "PrecompileTools", "Preferences", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"]
-git-tree-sha1 = "4891b745bd621f88aac661f2504d014931b443ba"
+git-tree-sha1 = "a5bc9c06e28108e04de0485273f0b5933cec66ed"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "3.10.0"
+version = "3.12.3"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveBandedMatricesExt = "BandedMatrices"
@@ -695,6 +798,12 @@ deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
 version = "0.8.1+2"
 
+[[deps.OpenSSL_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "3da7367955dcc5c54c1ba4d402ccdc09a1a3e046"
+uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
+version = "3.0.13+1"
+
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
@@ -706,6 +815,11 @@ git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.6.3"
 
+[[deps.PCRE2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
+version = "10.42.0+1"
+
 [[deps.PackageExtensionCompat]]
 git-tree-sha1 = "fb28e33b8a95c4cee25ce296c817d89cc2e53518"
 uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930"
@@ -719,10 +833,10 @@ uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 version = "0.12.3"
 
 [[deps.Parsers]]
-deps = ["Dates", "SnoopPrecompile"]
-git-tree-sha1 = "cceb0257b662528ecdf0b4b4302eb00e767b38e7"
+deps = ["Dates", "PrecompileTools", "UUIDs"]
+git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.5.0"
+version = "2.8.1"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -785,9 +899,9 @@ version = "1.3.4"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "758bc86b90e9fee2edc4af2a750b0d3f2d5c02c5"
+git-tree-sha1 = "d0f8d22294f932efb1617d669aff73a5c97d38ff"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.19.0"
+version = "3.20.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
@@ -818,6 +932,12 @@ git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
 version = "1.2.2"
 
+[[deps.RegistryInstances]]
+deps = ["LazilyInitializedFields", "Pkg", "TOML", "Tar"]
+git-tree-sha1 = "ffd19052caf598b8653b99404058fce14828be51"
+uuid = "2792f1a3-b283-48e8-9a74-f99dce5104f3"
+version = "0.1.0"
+
 [[deps.Requires]]
 deps = ["UUIDs"]
 git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
@@ -847,9 +967,9 @@ version = "0.6.42"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
-git-tree-sha1 = "265f1a7a804d8093fa0b17e33e45373a77e56ca5"
+git-tree-sha1 = "9f59654e2a85017ee27b0f59c7fac5a57aa10ced"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.38.0"
+version = "2.39.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -895,24 +1015,20 @@ deps = ["Distributed", "Mmap", "Random", "Serialization"]
 uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 
 [[deps.SimpleNonlinearSolve]]
-deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "DiffResults", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"]
-git-tree-sha1 = "c020028bb22a2f23cbd88cb92cf47cbb8c98513f"
+deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "DiffResults", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"]
+git-tree-sha1 = "f3c50acd5145f2c6ee85343ce6f433dd2caab41e"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-version = "1.8.0"
+version = "1.9.0"
 
     [deps.SimpleNonlinearSolve.extensions]
     SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore"
-    SimpleNonlinearSolvePolyesterForwardDiffExt = "PolyesterForwardDiff"
     SimpleNonlinearSolveReverseDiffExt = "ReverseDiff"
-    SimpleNonlinearSolveStaticArraysExt = "StaticArrays"
     SimpleNonlinearSolveTrackerExt = "Tracker"
     SimpleNonlinearSolveZygoteExt = "Zygote"
 
     [deps.SimpleNonlinearSolve.weakdeps]
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
     Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
@@ -922,11 +1038,6 @@ git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
 uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 version = "0.9.4"
 
-[[deps.SnoopPrecompile]]
-git-tree-sha1 = "f604441450a3c0569830946e5b33b78c928e1a85"
-uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
-version = "1.0.1"
-
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
@@ -937,9 +1048,9 @@ version = "1.10.0"
 
 [[deps.SparseDiffTools]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "PackageExtensionCompat", "Random", "Reexport", "SciMLOperators", "Setfield", "SparseArrays", "StaticArrayInterface", "StaticArrays", "Tricks", "UnPack", "VertexSafeGraphs"]
-git-tree-sha1 = "cce98ad7c896e52bb0eded174f02fc2a29c38477"
+git-tree-sha1 = "469f51f8c4741ce944be2c0b65423b518b1405b0"
 uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
-version = "2.18.0"
+version = "2.19.0"
 
     [deps.SparseDiffTools.extensions]
     SparseDiffToolsEnzymeExt = "Enzyme"
@@ -955,6 +1066,12 @@ version = "2.18.0"
     Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
     Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
+[[deps.SparseMatrixColorings]]
+deps = ["ADTypes", "Compat", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "d4adedbcc8776c567e0e22ef19f13cf10cb6ecaa"
+uuid = "0a514795-09f3-496d-8182-132a7b665d35"
+version = "0.3.2"
+
 [[deps.Sparspak]]
 deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"]
 git-tree-sha1 = "342cf4b449c299d8d1ceaf00b7a49f4fbc7940e7"
@@ -966,12 +1083,16 @@ deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_j
 git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.4.0"
+weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
-    [deps.SpecialFunctions.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+[[deps.StackViews]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
+uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
+version = "0.1.1"
 
 [[deps.Static]]
 deps = ["IfElse"]
@@ -995,15 +1116,12 @@ deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
 git-tree-sha1 = "9ae599cd7529cfce7fea36cf00a62cfc56f0f37c"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
 version = "1.9.4"
+weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
     StaticArraysStatisticsExt = "Statistics"
 
-    [deps.StaticArrays.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -1031,9 +1149,9 @@ version = "7.2.1+1"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "b479c7a16803f08779ac5b7f9844a42621baeeda"
+git-tree-sha1 = "a5f6f138b740c9d93d76f0feddd3092e6ef002b7"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.21"
+version = "0.3.22"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -1073,6 +1191,15 @@ git-tree-sha1 = "5a13ae8a41237cff5ecf34f73eb1b8f42fff6531"
 uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 version = "0.5.24"
 
+[[deps.TranscodingStreams]]
+git-tree-sha1 = "5d54d076465da49d6746c647022f3b3674e64156"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.10.8"
+weakdeps = ["Random", "Test"]
+
+    [deps.TranscodingStreams.extensions]
+    TestExt = ["Test", "Random"]
+
 [[deps.TriangularSolve]]
 deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"]
 git-tree-sha1 = "66c68a20907800c0b7c04ff8a6164115e8747de2"
@@ -1104,9 +1231,9 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [[deps.VectorizationBase]]
 deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
-git-tree-sha1 = "6129a4faf6242e7c3581116fbe3270f3ab17c90d"
+git-tree-sha1 = "e863582a41c5731f51fd050563ae91eb33cf09be"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.21.67"
+version = "0.21.68"
 
 [[deps.VertexSafeGraphs]]
 deps = ["Graphs"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 86d55f6..8064d77 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 NICE = "170674ed-ea2e-4c97-a263-9f057371600a"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
diff --git a/docs/make.jl b/docs/make.jl
index b256246..9dd3216 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -9,6 +9,7 @@ makedocs(;
     repo="https://github.com/theochem/NICE.jl/blob/{commit}{path}#{line}",
     sitename="NICE.jl",
     format=Documenter.HTML(;
+        repolink="https://github.com/theochem/NICE.jl/",
         prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://nice.qcdevs.org/",
         edit_link="main",
diff --git a/docs/src/index.md b/docs/src/index.md
index dd9b3d9..ccc661a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -4,7 +4,7 @@ CurrentModule = NICE
 
 # NICE.jl
 
-Documentation for [NICE.jl](nice.qcdevs.org).
+Documentation for NICE.jl.
 
 ```@contents
 Pages = ["api.md"]
diff --git a/src/Exact.jl b/src/Exact.jl
index a020989..1333bdd 100644
--- a/src/Exact.jl
+++ b/src/Exact.jl
@@ -1,11 +1,25 @@
 using LinearAlgebra
 
+using Reexport
+
 using SciMLBase
 using ForwardDiff
 using NonlinearSolve
 
+using Evolutionary
+@reexport using Evolutionary: ES, CMAES, GA, DE, NSGA2, TreeGP
+
+export solve!
+export evolutionary_solve!
+
+struct ConvergenceError <: Exception
+    msg::String
+end
+
+Base.showerror(io::IO, e::ConvergenceError) = print(io, e.msg)
+
 """
-    solve(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
+    solve!(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
 
 Solve the system deterministically.
 
@@ -21,58 +35,42 @@ K_\\text{eq}^{\\left(m\\right)} =
     \\prod_{n=1}^N {\\left( a_n^{\\left(0\\right)} +
     \\sum_{m=1}^M c_{mn} \\xi_m \\right)}^{c_{mn}}
 ```
-using user-provided ``K_\text{eq}`` values and the Newton's method + Trust
+using user-provided ``K_\\text{eq}`` values and the Newton's method + Trust
 Region method.
 """
-function solve(
+function solve!(
     rxn_system::ReactionSystem,
-    K_eqs::AbstractVector{Float64};
-    maxiters::Integer=1000,
-    abstol::Real=1.0e-9,
-    reltol::Real=0.0,
+    K_eqs::Vector{Float64};
+    maxiter::Int=1000,
+    abstol::Float64=1.0e-9,
+    reltol::Float64=0.0,
 )
-    # Define the objective function
-    function f_K_eqs(f, ξ, _)
-        for i in 1 : rxn_system.n_reaction
-            f_i1 = K_eqs[i]
-            f_i2 = 1.0
-            for j in 1 : rxn_system.n_species
-                # Reactants
-                if rxn_system.stoich[j, i] < 0
-                    f_i1 *= (
-                        rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
-                    ) ^ -rxn_system.stoich[j, i]
-                # Products
-                elseif rxn_system.stoich[j, i] > 0
-                    f_i2 *= (
-                        rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
-                    ) ^ rxn_system.stoich[j, i]
-                end
-            end
-            f[i] = f_i1 - f_i2
-        end
-        nothing
-    end
-    # Compute reaction extents
-    ξ = rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init)
-    # Create objective function and Jacobian (automatic, with autodiff)
+    # Define nonlinear problem (derives Jacobian via autodiff)
     problem = NonlinearSolve.NonlinearProblem(
-        f_K_eqs, ξ;
-        maxiters=maxiters, abstol=abstol, reltol=reltol,
+        # Objective function
+        (f, ξ, _) -> f_K_eqs(f, ξ, rxn_system, K_eqs),
+        # Initial guess for reaction extents ξ
+        rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init);
+        # Solver options
+        maxiters=maxiter,
+        abstol=abstol,
+        reltol=reltol,
     )
     # Run the nonlinear optimization
     solution = NonlinearSolve.solve(problem, NonlinearSolve.TrustRegion())
-    # If nonlinear optimization was successful, update `rxn_system.concs`
-    if SciMLBase.successful_retcode(solution)
-        rxn_system.concs .= rxn_system.concs_init
-        rxn_system.concs .+= rxn_system.stoich * solution.u
+    # If nonlinear optimization was successful, update `rxn_system.concs`;
+    # otherwise throw ConvergenceError exception
+    if !SciMLBase.successful_retcode(solution)
+        throw(ConvergenceError("Did not converge"))
     end
+    rxn_system.concs .= rxn_system.concs_init
+    rxn_system.concs .+= rxn_system.stoich * solution.u
     # Return the nonlinear solution
     solution
 end
 
 """
-    solve(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
+    solve!(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
 
 Solve the system deterministically.
 
@@ -94,13 +92,100 @@ Instead of specifying the ``K_\\text{eq}`` values directly, they are
 approximated here using the forward or reverse rate constants and a parameter
 ``\\phi``.
 """
-function solve(
+function solve!(
     rxn_system::ReactionSystem;
-    φ::Real=1.0,
+    φ::Float64=1.0,
     rate_consts::Symbol=:forward,
-    maxiters::Integer=1000,
-    abstol::Real=1.0e-9,
-    reltol::Real=0.0,
+    maxiter::Int=1000,
+    abstol::Float64=1.0e-9,
+    reltol::Float64=0.0,
+)
+    solve!(
+        rxn_system,
+        approximate_K_eqs(rxn_system, φ, rate_consts);
+        maxiter=maxiter,
+        abstol=abstol,
+        reltol=reltol,
+    )
+end
+
+"""
+    evolutionary_solve!(rxn_system, K_eqs, algorithm; abstol=Inf, reltol=Inf, maxiter=1000)
+
+Solve the system using an evolutionary algorithm.
+
+Solves the nonlinear problem described by [`solve!`](@ref) by using one of the
+evolutionary algorithms from
+[Evolutionary.jl](https://wildart.github.io/Evolutionary.jl/stable/).
+These algorithm types are reexported by NICE.jl.
+"""
+function evolutionary_solve!(
+    rxn_system::ReactionSystem,
+    K_eqs::Vector{Float64},
+    algorithm::TOpt;
+    abstol::Float64=Inf,
+    reltol::Float64=Inf,
+    maxiter::Int=1000,
+) where {TOpt <: Evolutionary.AbstractOptimizer}
+    f = similar(K_eqs)
+    solution = Evolutionary.optimize(
+        # Objective function
+        ξ -> f_K_eqs_cma(f, ξ, rxn_system, K_eqs),
+        # Initial guess for reaction extents ξ
+        rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init),
+        # Evolutionary algorithm
+        algorithm,
+        # Solver options
+        Evolutionary.Options(iterations=maxiter, abstol=abstol, reltol=reltol),
+    )
+    # If nonlinear optimization was successful, update `rxn_system.concs`
+    # otherwise throw ConvergenceError exception
+    if !solution.converged
+        throw(ConvergenceError("Did not converge"))
+    end
+    rxn_system.concs .= rxn_system.concs_init
+    rxn_system.concs .+= rxn_system.stoich * solution.minimizer
+    # Return the nonlinear solution
+    solution
+end
+
+"""
+    evolutionary_solve!(rxn_system, K_eqs, algorithm; abstol=Inf, reltol=Inf, maxiter=1000)
+
+Solve the system using an evolutionary algorithm.
+
+Solves the nonlinear problem described by [`solve!`](@ref) by using one of the
+evolutionary algorithms from
+[Evolutionary.jl](https://wildart.github.io/Evolutionary.jl/stable/).
+These algorithm types are reexported by NICE.jl.
+
+Instead of specifying the ``K_\\text{eq}`` values directly, they are
+approximated here using the forward or reverse rate constants and a parameter
+``\\phi``.
+"""
+function evolutionary_solve!(
+    rxn_system::ReactionSystem,
+    algorithm::TOpt;
+    φ::Float64=1.0,
+    rate_consts::Symbol=:forward,
+    abstol::Float64=Inf,
+    reltol::Float64=Inf,
+    maxiter::Int=1000,
+) where {TOpt <: Evolutionary.AbstractOptimizer}
+    evolutionary_solve!(
+        rxn_system,
+        approximate_K_eqs(rxn_system, φ, rate_consts),
+        algorithm;
+        abstol=abstol,
+        reltol=reltol,
+        maxiter=maxiter,
+    )
+end
+
+function approximate_K_eqs(
+    rxn_system::ReactionSystem,
+    φ::Float64,
+    rate_consts::Symbol,
 )
     if rate_consts === :forward
         # Evaluate K_eqs from forward rate constants and φ_fwd
@@ -115,5 +200,32 @@ function solve(
         throw(ArgumentError("invalid rate_consts value $rate_consts \
                              (must be :forward or :reverse)"))
     end
-    solve(rxn_system, K_eqs; maxiters=maxiters, abstol=abstol, reltol=reltol)
+    K_eqs
+end
+
+function f_K_eqs(f, ξ, rxn_system, K_eqs)::Nothing
+    for i in 1 : rxn_system.n_reaction
+        f_i1 = K_eqs[i]
+        f_i2 = 1.0
+        for j in 1 : rxn_system.n_species
+            # Reactants
+            if rxn_system.stoich[j, i] < 0
+                f_i1 *= (
+                    rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
+                ) ^ -rxn_system.stoich[j, i]
+            # Products
+            elseif rxn_system.stoich[j, i] > 0
+                f_i2 *= (
+                    rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
+                ) ^ rxn_system.stoich[j, i]
+            end
+        end
+        f[i] = f_i1 - f_i2
+    end
+    nothing
+end
+
+function f_K_eqs_cma(f::Vector{Float64}, ξ::Vector{Float64}, rxn_system::ReactionSystem, K_eqs::Vector{Float64})::Float64
+    f_K_eqs(f, ξ, rxn_system, K_eqs)
+    √(f · f)
 end
diff --git a/src/NEKMC.jl b/src/NEKMC.jl
index c9d4d67..b0c5b5c 100644
--- a/src/NEKMC.jl
+++ b/src/NEKMC.jl
@@ -1,7 +1,9 @@
 using LinearAlgebra
 
+export simulate!
+
 """
-    simulate(rxn_system;
+    simulate!(rxn_system;
              n_iter=Int(1e+8), chunk_iter=Int(1e+4),
              ε=1.0e-4, ε_scale=1.0, ε_concs=0.0,
              tol_t=Inf, tol_ε=0.0, tol_concs=0.0)
@@ -9,16 +11,16 @@ using LinearAlgebra
 Run a *N*et-*E*vent *K*inetic *M*onte *C*arlo (NEKMC) simulation to find the
 equilibrium concentrations of the reaction system.
 """
-function simulate(
+function simulate!(
     rxn_system::ReactionSystem;
-    n_iter::Integer=Int(1e+8),
-    chunk_iter::Integer=Int(1e+4),
-    ε::Real=1.0e-4,
-    ε_scale::Real=1.0,
-    ε_concs::Real=0.0,
-    tol_t::Real=Inf,
-    tol_ε::Real=0.0,
-    tol_concs::Real=0.0,
+    n_iter::Int=Int(1e+8),
+    chunk_iter::Int=Int(1e+4),
+    ε::Float64=1.0e-4,
+    ε_scale::Float64=1.0,
+    ε_concs::Float64=0.0,
+    tol_t::Float64=Inf,
+    tol_ε::Float64=0.0,
+    tol_concs::Float64=0.0,
 )
     # Allocate probability vector, Δconcs vector, Δt
     pvec = zeros(Float64, rxn_system.n_reaction)
@@ -82,7 +84,7 @@ end
 
 function select_reaction(
     rxn_system::ReactionSystem,
-    pvec::AbstractVector{Float64},
+    pvec::Vector{Float64},
 )
     # Update probability vector
     p = 0.
@@ -103,8 +105,8 @@ end
 
 function do_reaction(
     rxn_system::ReactionSystem,
-    i_rxn::Integer,
-    ε::Real,
+    i_rxn::Int,
+    ε::Float64,
 )
     rate = rxn_system.net_rates[i_rxn]
     if rate >= 0
diff --git a/src/NICE.jl b/src/NICE.jl
index 3c80c28..858e3a6 100644
--- a/src/NICE.jl
+++ b/src/NICE.jl
@@ -1,16 +1,13 @@
 """
 NICE (N-species Ice Table) module -- provides simultaneous equilibria solvers:
-- [`simulate`](@ref) Net-event kinetic Monte Carlo solver
-- [`solve`](@ref) Nonlinear equations solver (Newton trust region method)
+- [`simulate!`](@ref) Net-event kinetic Monte Carlo solver
+- [`solve!`](@ref) Exact nonlinear equations solver (Newton trust region method)
+- [`evolutionary_solve!`](@ref) Evolutionary nonlinear equations solver (Newton trust region method)
 
 The [`ReactionSystem`](@ref) type is provided as an interface for these solvers.
 """
 module NICE
 
-export ReactionSystem
-export simulate
-export solve
-
 include("ReactionSystem.jl")
 include("NEKMC.jl")
 include("Exact.jl")
diff --git a/src/ReactionSystem.jl b/src/ReactionSystem.jl
index 6294ba8..16ce029 100644
--- a/src/ReactionSystem.jl
+++ b/src/ReactionSystem.jl
@@ -1,7 +1,9 @@
+export ReactionSystem
+
 """
 Definition of the simultaneous equilibria reaction system. This type is the
-input to the solver functions [`simulate`](@ref) and [`solve`](@ref), and is
-updated in-place.
+input to the solver functions [`simulate!`](@ref), [`solve!`](@ref) and
+[`evolutionary_solve!`](@ref), and is updated in-place.
 
 The reaction system consists of
 - ``M`` Reactions involving ``N`` species;
@@ -50,10 +52,10 @@ Instantiate a [`ReactionSystem`](@ref) instance using the forward and reverse
 rate constants for each reaction.
 """
 function ReactionSystem(
-    stoich::AbstractMatrix{Float64},
-    concs::AbstractVector{Float64},
-    rev_rate_consts::AbstractVector{Float64},
-    fwd_rate_consts::AbstractVector{Float64},
+    stoich::Matrix{Float64},
+    concs::Vector{Float64},
+    rev_rate_consts::Vector{Float64},
+    fwd_rate_consts::Vector{Float64},
 )
     n_species = size(stoich, 1)
     n_reaction = size(stoich, 2)
@@ -84,10 +86,10 @@ k_f &= \\phi - k_r
 ```
 """
 function ReactionSystem(
-    stoich::AbstractMatrix{Float64},
-    concs::AbstractVector{Float64},
-    K_eqs::AbstractVector{Float64};
-    φ::Real=1.0,
+    stoich::Matrix{Float64},
+    concs::Vector{Float64},
+    K_eqs::Vector{Float64};
+    φ::Float64=1.0,
 )
     n_species = size(stoich, 1)
     n_reaction = size(stoich, 2)
diff --git a/test/runtests.jl b/test/runtests.jl
index 6a9c74b..3393d93 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,41 +1,82 @@
 using NICE
 using Test
 
-@testset "System 1" begin
-
-    test_concs = [0.9261879203, 0.9623865752, 0.0926187920]
-    stoich = [-0.5 -0.5;
-               1.0 -1.0;
-               0.0  1.0]
-    keq_vals = [1.0, 0.1]
-
-    concs = [1.0, 0.2, 0.4]
-    rxn_system = ReactionSystem(stoich, concs, keq_vals; φ=1.0)
-    simulate(rxn_system; n_iter=Int(5e+8), chunk_iter=Int(1e+3), ε=1e-3, ε_scale=0.5, ε_concs=2.0, tol_ε=0.)
-    @test isapprox(rxn_system.concs, test_concs; atol=1.0e-5)
-
-    solve(rxn_system, keq_vals)
-    @test isapprox(rxn_system.concs, test_concs; atol=1.0e-6)
-
+struct TestSystem
+    stoich::Matrix{Float64}
+    K_eqs::Vector{Float64}
+    concs::Vector{Float64}
+    concs_eq::Vector{Float64}
 end
 
-@testset "System 2" begin
-
-    test_concs = [9.0e-01, 1.10212630e-08, 1.00002433e-10, 9.98999891e-02, 0.0, 9.99999e-05]
-    stoich = [-1.0  0.0  0.0 -1.0;
-              -1.0 -1.0  0.0  0.0;
-               0.0 -1.0 -1.0  0.0;
-               1.0  0.0 -1.0  0.0;
-               0.0  1.0  0.0 -1.0;
-               0.0  0.0  1.0  1.0]
-    keq_vals = [1.0e7, 1.0e9, 1.0e7, 1.0e9]
+system1 = TestSystem(
+    # stoich
+    [-0.5 -0.5;
+      1.0 -1.0;
+      0.0  1.0],
+    # K_eqs
+    [1.0, 0.1],
+    # concs
+    [1.0, 0.2, 0.4],
+    # concs_eq
+    [0.9261879203, 0.9623865752, 0.0926187920],
+)
 
-    concs = [1.0, 0.1, 1e-4, 0.0, 0.0, 0.0]
-    rxn_system = ReactionSystem(stoich, concs, keq_vals; φ=1.0)
-    simulate(rxn_system; n_iter=Int(5e+8), chunk_iter=Int(1e+3), ε=1e-3, ε_scale=0.5, ε_concs=2.0, tol_ε=0.)
-    @test isapprox(rxn_system.concs, test_concs; atol=1.0e-5)
+system2 = TestSystem(
+    # stoich
+    [-1.0  0.0  0.0 -1.0;
+     -1.0 -1.0  0.0  0.0;
+      0.0 -1.0 -1.0  0.0;
+      1.0  0.0 -1.0  0.0;
+      0.0  1.0  0.0 -1.0;
+      0.0  0.0  1.0  1.0],
+    # K_eqs
+    [1.0e7, 1.0e9, 1.0e7, 1.0e9],
+    # concs
+    [1.0, 0.1, 1e-4, 0.0, 0.0, 0.0],
+    # concs_eq
+    [9.0e-01, 1.10212630e-08, 1.00002433e-10, 9.98999891e-02, 0.0, 9.99999e-05],
+)
 
-    solve(rxn_system, keq_vals)
-    @test isapprox(rxn_system.concs, test_concs; atol=1.0e-6)
+@testset verbose=true "System 1" begin
+    # Set up reaction system
+    rxn_system1 = ReactionSystem(system1.stoich, copy(system1.concs), system1.K_eqs; φ=1.0)
+    # Run and test NEKMC
+    simulate!(
+        rxn_system1;
+        n_iter=Int(1e+9),
+        chunk_iter=Int(1e+3),
+        ε=1e-3,
+        ε_scale=0.5,
+        ε_concs=2.0,
+        tol_ε=0.0,
+    )
+    @test isapprox(rxn_system1.concs, system1.concs_eq; rtol=1.0e-4)
+    # Run and test exact nonlinear solver
+    rxn_system1_copy = deepcopy(rxn_system1)
+    solve!(rxn_system1_copy, system1.K_eqs; reltol=1.0e-6)
+    @test isapprox(rxn_system1_copy.concs, system1.concs_eq; rtol=1.0e-6)
+    # Run and test CMA solver
+    rxn_system1_copy = deepcopy(rxn_system1)
+    evolutionary_solve!(rxn_system1_copy, system1.K_eqs, CMAES(); reltol=1.0e-6)
+    @test isapprox(rxn_system1_copy.concs, system1.concs_eq; rtol=1.0e-6)
+end
 
+@testset verbose=true "System 2" begin
+    # Set up reaction system
+    rxn_system2 = ReactionSystem(system2.stoich, copy(system2.concs), system2.K_eqs; φ=1.0)
+    # Run and test NEKMC
+    simulate!(
+        rxn_system2;
+        n_iter=Int(1e+9),
+        chunk_iter=Int(1e+2),
+        ε=1e-5,
+        ε_scale=0.5,
+        ε_concs=2.0,
+        # tol_ε=0.0,
+    )
+    @test isapprox(rxn_system2.concs, system2.concs_eq; rtol=1.0e-4)
+    # Run and test CMA solver
+    rxn_system2_copy = deepcopy(rxn_system2)
+    evolutionary_solve!(rxn_system2_copy, system2.K_eqs, CMAES(); maxiter=1000000, reltol=1.0e-6)
+    @test isapprox(rxn_system2_copy.concs, system2.concs_eq; rtol=1.0e-6)
 end