diff --git a/benchmarks/NonStiffODE/LotkaVolterra_wpd.jmd b/benchmarks/NonStiffODE/LotkaVolterra_wpd.jmd
index cd06ec4d6..c50de34be 100644
--- a/benchmarks/NonStiffODE/LotkaVolterra_wpd.jmd
+++ b/benchmarks/NonStiffODE/LotkaVolterra_wpd.jmd
@@ -9,7 +9,7 @@ The purpose of this problem is to test the performance on easy problems. Since i
 
 ```julia
 using OrdinaryDiffEq, ParameterizedFunctions, ODEInterfaceDiffEq, LSODA,
-      Sundials, DiffEqDevTools, StaticArrays, OrdinaryDiffEqSIMDRK
+      Sundials, DiffEqDevTools, StaticArrays, OrdinaryDiffEqSIMDRK, IRKGaussLegendre
 
 f = @ode_def LotkaVolterra begin
     dx = a*x - b*x*y
@@ -42,6 +42,8 @@ setups = [Dict(:alg=>DP5())
           Dict(:alg=>dopri5())
           Dict(:alg=>Tsit5())
           Dict(:alg=>Vern6())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Tsit5(), :prob_choice => 2)
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>MER5v2(), :prob_choice => 2)
@@ -63,9 +65,12 @@ setups = [Dict(:alg=>DP5())
           #Dict(:alg=>ode45())
           Dict(:alg=>Tsit5())
           Dict(:alg=>Vern6())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Tsit5(), :prob_choice => 2)
           Dict(:alg=>Vern6(), :prob_choice => 2)]
-wp = WorkPrecisionSet(probs, abstols, reltols, setups; appxsol = test_sol, maxiters = 10000,
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups; appxsol = test_sol, maxiters = 10000,
     error_estimate = :L2, dense_errors = true, numruns = 100)
 plot(wp)
 ```
@@ -82,6 +87,8 @@ setups = [Dict(:alg=>DP8())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>Vern7(), :prob_choice => 2)
           Dict(:alg=>Vern8(), :prob_choice => 2)
@@ -98,6 +105,8 @@ setups = [Dict(:alg=>odex())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>Vern7(), :prob_choice => 2)
           Dict(:alg=>Vern8(), :prob_choice => 2)
@@ -119,6 +128,8 @@ setups = [Dict(:alg=>DP8())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>Vern7(), :prob_choice => 2)
           Dict(:alg=>Vern8(), :prob_choice => 2)
@@ -141,6 +152,8 @@ reltols = 1.0 ./ 10.0 .^ (8:13)
 setups = [Dict(:alg=>Tsit5())
           Dict(:alg=>Vern9())
           Dict(:alg=>VCABM())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern9(), :prob_choice => 2)
           Dict(:alg=>VCABM(), :prob_choice => 2)
           Dict(:alg=>AitkenNeville(min_order = 1, max_order = 9, init_order = 4, threading = true))
@@ -148,7 +161,8 @@ setups = [Dict(:alg=>Tsit5())
               min_order = 1, max_order = 9, init_order = 4, threading = true))
           Dict(:alg=>ExtrapolationMidpointHairerWanner(
               min_order = 2, max_order = 11, init_order = 4, threading = true))]
-solnames = ["Tsit5", "Vern9", "VCABM", "Vern9 Static", "VCABM Static",
+solnames = ["Tsit5", "Vern9", "VCABM", "IRKGL16", "IRKGL16 SIMD",
+    "Vern9 Static", "VCABM Static",
     "AitkenNeville", "Midpoint Deuflhard", "Midpoint Hairer Wanner"]
 wp = WorkPrecisionSet(
     probs, abstols, reltols, setups; appxsol = test_sol, names = solnames,
diff --git a/benchmarks/NonStiffODE/Manifest.toml b/benchmarks/NonStiffODE/Manifest.toml
index 22597045c..e0085b9df 100644
--- a/benchmarks/NonStiffODE/Manifest.toml
+++ b/benchmarks/NonStiffODE/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.10"
 manifest_format = "2.0"
-project_hash = "da2771e14f14e103639884677f2c62dad75406b5"
+project_hash = "5914eb28623bb0e6c136445a1023af008caad763"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "60665b326b75db6517939d0e1875850bc4a54368"
@@ -73,9 +73,9 @@ version = "0.4.0"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "9606d7832795cbef89e06a550475be300364a8aa"
+git-tree-sha1 = "dbd8c3bbbdbb5c2778f85f4422c39960eac65a42"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.19.0"
+version = "7.20.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -85,6 +85,7 @@ version = "7.19.0"
     ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
     ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceMetalExt = "Metal"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceSparseArraysExt = "SparseArrays"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
@@ -98,6 +99,7 @@ version = "7.19.0"
     ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -137,9 +139,9 @@ version = "0.1.6"
 
 [[deps.BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "84a4360c718e7473fec971ae27f409a2f24befc8"
+git-tree-sha1 = "79e651aa489a7879107d66e3d1948e9aa1b4055e"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "1.7.1"
+version = "1.7.2"
 
     [deps.BlockArrays.extensions]
     BlockArraysAdaptExt = "Adapt"
@@ -240,9 +242,9 @@ uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.13.1"
 
 [[deps.Combinatorics]]
-git-tree-sha1 = "8010b6bb3388abe68d95743dcbea77650bb2eddf"
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
 uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-version = "1.0.3"
+version = "1.0.2"
 
 [[deps.CommonMark]]
 deps = ["PrecompileTools"]
@@ -368,14 +370,13 @@ version = "1.9.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "529087a10b502502201c7bf4e8769e40d7c15a17"
+git-tree-sha1 = "c969d82ad56c5272cf670e0c17f75420f1668ac7"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.184.0"
+version = "6.188.0"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
     DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
-    DiffEqBaseDistributionsExt = "Distributions"
     DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
     DiffEqBaseForwardDiffExt = ["ForwardDiff"]
     DiffEqBaseGTPSAExt = "GTPSA"
@@ -588,9 +589,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.5"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "8272a687bca7b5c601c0c24fc0c71bff10aafdfd"
+git-tree-sha1 = "787f5e2efb33af12cc3fde313c1b0195a3f08eba"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.12"
+version = "0.8.13"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -749,9 +750,9 @@ version = "1.3.7"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "a2df1b776752e3f344e5116c06d75a10436ab853"
+git-tree-sha1 = "ce15956960057e9ff7f1f535400ffa14c92429a4"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.38"
+version = "1.1.0"
 weakdeps = ["StaticArrays"]
 
     [deps.ForwardDiff.extensions]
@@ -826,6 +827,12 @@ git-tree-sha1 = "45288942190db7c5f760f59c04495064eedf9340"
 uuid = "b0724c58-0f36-5564-988d-3bb0596ebc4a"
 version = "0.22.4+0"
 
+[[deps.Ghostscript_jll]]
+deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "38044a04637976140074d0b0621c1edf0eb531fd"
+uuid = "61579ee1-b43e-5ca0-a5da-69d92c66a64b"
+version = "9.55.1+0"
+
 [[deps.Git]]
 deps = ["Git_LFS_jll", "Git_jll", "JLLWrappers", "OpenSSH_jll"]
 git-tree-sha1 = "824a1890086880696fc908fe12a17bcf61738bd8"
@@ -840,9 +847,9 @@ version = "3.7.0+0"
 
 [[deps.Git_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "cd06e503111a7c5ef1d4a339de6ccf5bd7437b32"
+git-tree-sha1 = "e2aef26f7d273f1e5b1daba56837c47b49b4388f"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
-version = "2.51.0+0"
+version = "2.51.1+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
@@ -904,9 +911,15 @@ version = "0.3.28"
 
 [[deps.IJulia]]
 deps = ["Base64", "Conda", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Pkg", "PrecompileTools", "Printf", "REPL", "Random", "SHA", "Sockets", "UUIDs", "ZMQ"]
-git-tree-sha1 = "66b1f8d8b870814a079da381dd20416c53e73839"
+git-tree-sha1 = "00a9cbc52b819846e89ef5b221227c33b110d8e0"
 uuid = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
-version = "1.30.1"
+version = "1.30.3"
+
+[[deps.IRKGaussLegendre]]
+deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "Parameters", "RecursiveArrayTools", "Reexport", "SIMD", "SciMLBase"]
+git-tree-sha1 = "5bcafc076844aa7055506775320332ccdfc36498"
+uuid = "58bc7355-f626-4c51-96f2-1f8a038f95a2"
+version = "0.2.11"
 
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
@@ -996,9 +1009,9 @@ version = "0.2.1"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e95866623950267c1e4878846f848d94810de475"
+git-tree-sha1 = "4255f0032eafd6451d707a51d5f0248b8a165e4d"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "3.1.2+0"
+version = "3.1.3+0"
 
 [[deps.JuliaFormatter]]
 deps = ["CommonMark", "Glob", "JuliaSyntax", "PrecompileTools", "TOML"]
@@ -1073,10 +1086,10 @@ uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 version = "1.4.0"
 
 [[deps.Latexify]]
-deps = ["Format", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Requires"]
-git-tree-sha1 = "52e1296ebbde0db845b356abbbe67fb82a0a116c"
+deps = ["Format", "Ghostscript_jll", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Requires"]
+git-tree-sha1 = "44f93c47f9cd6c7e431f2f2091fcba8f01cd7e8f"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-version = "0.16.9"
+version = "0.16.10"
 
     [deps.Latexify.extensions]
     DataFramesExt = "DataFrames"
@@ -1203,9 +1216,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "OpenBLAS_jll", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "StaticArraysCore", "UnPack"]
-git-tree-sha1 = "0f1a02cea457a2e26b67e105aa7ee549419c2550"
+git-tree-sha1 = "3636b906db79eb1c5dbd97c234cbd8e7a08721dd"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.37.0"
+version = "3.39.2"
 
     [deps.LinearSolve.extensions]
     LinearSolveAMDGPUExt = "AMDGPU"
@@ -1610,9 +1623,9 @@ version = "1.10.1"
 
 [[deps.OrdinaryDiffEqCore]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "e579c9a4f9102e82da3d97c349a74d6bc11cf8dc"
+git-tree-sha1 = "8d50e4c5cd9de05602e2fa38161e896454a13c20"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-version = "1.30.0"
+version = "1.33.0"
 
     [deps.OrdinaryDiffEqCore.extensions]
     OrdinaryDiffEqCoreEnzymeCoreExt = "EnzymeCore"
@@ -1629,10 +1642,14 @@ uuid = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
 version = "1.8.0"
 
 [[deps.OrdinaryDiffEqDifferentiation]]
-deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"]
-git-tree-sha1 = "4c270747152db513dd80bfd5f2f9df48befff28a"
+deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"]
+git-tree-sha1 = "1fff7356ccb79ca1b1ff3276135fda68d4d661be"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
-version = "1.14.0"
+version = "1.16.0"
+weakdeps = ["SparseArrays"]
+
+    [deps.OrdinaryDiffEqDifferentiation.extensions]
+    OrdinaryDiffEqDifferentiationSparseArraysExt = "SparseArrays"
 
 [[deps.OrdinaryDiffEqExplicitRK]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
@@ -1648,9 +1665,9 @@ version = "1.8.0"
 
 [[deps.OrdinaryDiffEqExtrapolation]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "FastPower", "LinearSolve", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "ee2cba2533e9faf71b09a319a910d4886931e7a6"
+git-tree-sha1 = "9e1b11cf448a2c1bca640103c1c848a20aa2f967"
 uuid = "becaefa8-8ca2-5cf9-886d-c06f3d2bd2c4"
-version = "1.8.0"
+version = "1.9.0"
 
 [[deps.OrdinaryDiffEqFIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "FastGaussQuadrature", "FastPower", "LinearAlgebra", "LinearSolve", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators"]
@@ -1696,9 +1713,9 @@ version = "1.6.0"
 
 [[deps.OrdinaryDiffEqLowStorageRK]]
 deps = ["Adapt", "DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "StaticArrays"]
-git-tree-sha1 = "9291cdfd2e8c91e900c48d71d76618de47daeede"
+git-tree-sha1 = "e6bd0a7fb6643a57b06a90415608a81aaf7bd772"
 uuid = "b0944070-b475-4768-8dec-fb6eb410534d"
-version = "1.6.0"
+version = "1.7.0"
 
 [[deps.OrdinaryDiffEqNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "FastClosures", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "PreallocationTools", "RecursiveArrayTools", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "StaticArrays"]
@@ -1738,9 +1755,9 @@ version = "1.5.0"
 
 [[deps.OrdinaryDiffEqRosenbrock]]
 deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "a06c1263d71ea42a1881b4d49c8a087035d4a3ff"
+git-tree-sha1 = "d0b4e34792fb64c3815fc79ad3adc300b1e35588"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "1.16.1"
+version = "1.17.0"
 
 [[deps.OrdinaryDiffEqSDIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
@@ -1756,9 +1773,9 @@ version = "1.1.0"
 
 [[deps.OrdinaryDiffEqSSPRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "StaticArrays"]
-git-tree-sha1 = "af955c61407631d281dd4c2e8331cdfea1af49be"
+git-tree-sha1 = "3bce87977264916bd92455754ab336faec68bf8a"
 uuid = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
-version = "1.6.0"
+version = "1.7.0"
 
 [[deps.OrdinaryDiffEqStabilizedIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqStabilizedRK", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays"]
@@ -1894,9 +1911,9 @@ version = "0.2.4"
 
 [[deps.PreallocationTools]]
 deps = ["Adapt", "ArrayInterface", "PrecompileTools"]
-git-tree-sha1 = "9b4ee15d1fc68654031964a7af0c914c898e35a7"
+git-tree-sha1 = "c05b4c6325262152483a1ecb6c69846d2e01727b"
 uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-version = "0.4.33"
+version = "0.4.34"
 
     [deps.PreallocationTools.extensions]
     PreallocationToolsForwardDiffExt = "ForwardDiff"
@@ -2096,6 +2113,12 @@ version = "1.4.0"
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
+[[deps.SIMD]]
+deps = ["PrecompileTools"]
+git-tree-sha1 = "fea870727142270bdf7624ad675901a1ee3b4c87"
+uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
+version = "3.7.1"
+
 [[deps.SIMDTypes]]
 git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
 uuid = "94e857df-77ce-4151-89e5-788b33177be4"
@@ -2109,9 +2132,9 @@ version = "0.6.43"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "33786c33502a6652fba239d3062ccc5e4cdd30c4"
+git-tree-sha1 = "3796862d75d9c5c5542bde5bf52ea932320c88f4"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.114.0"
+version = "2.116.3"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2119,6 +2142,7 @@ version = "2.114.0"
     SciMLBaseForwardDiffExt = "ForwardDiff"
     SciMLBaseMLStyleExt = "MLStyle"
     SciMLBaseMakieExt = "Makie"
+    SciMLBaseMeasurementsExt = "Measurements"
     SciMLBaseMonteCarloMeasurementsExt = "MonteCarloMeasurements"
     SciMLBaseMooncakeExt = "Mooncake"
     SciMLBasePartialFunctionsExt = "PartialFunctions"
@@ -2155,15 +2179,15 @@ version = "0.1.3"
 
 [[deps.SciMLJacobianOperators]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DifferentiationInterface", "FastClosures", "LinearAlgebra", "SciMLBase", "SciMLOperators"]
-git-tree-sha1 = "3414071e3458f3065de7fa5aed55283b236b4907"
+git-tree-sha1 = "a273b291c90909ba6fe08402dd68e09aae423008"
 uuid = "19f34311-ddf3-4b8b-af20-060888a46c0e"
-version = "0.1.8"
+version = "0.1.11"
 
 [[deps.SciMLOperators]]
 deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "MacroTools"]
-git-tree-sha1 = "78ac1b947205b07973321f67f17df8fbe6154ac9"
+git-tree-sha1 = "753026a73e1bee963239016a162db8dc1935a153"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.6.0"
+version = "1.7.1"
 weakdeps = ["SparseArrays", "StaticArraysCore"]
 
     [deps.SciMLOperators.extensions]
@@ -2305,9 +2329,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"]
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "cbea8a6bd7bed51b1619658dec70035e07b8502f"
+git-tree-sha1 = "b8693004b385c842357406e3af647701fe783f98"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.14"
+version = "1.9.15"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -2421,9 +2445,9 @@ version = "0.2.2"
 
 [[deps.SymbolicUtils]]
 deps = ["AbstractTrees", "ArrayInterface", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "ExproniconLite", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "TaskLocalValues", "TermInterface", "TimerOutputs", "Unityper"]
-git-tree-sha1 = "8c103c491ccf3e2b4284635c24b5de768adc6be8"
+git-tree-sha1 = "a85b4262a55dbd1af39bb6facf621d79ca6a322d"
 uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
-version = "3.31.0"
+version = "3.32.0"
 
     [deps.SymbolicUtils.extensions]
     SymbolicUtilsLabelledArraysExt = "LabelledArrays"
@@ -2435,9 +2459,9 @@ version = "3.31.0"
 
 [[deps.Symbolics]]
 deps = ["ADTypes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "LaTeXStrings", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "OffsetArrays", "PrecompileTools", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "SciMLPublic", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"]
-git-tree-sha1 = "5f2f0188931997573d2a8e24f43bac699da95d8f"
+git-tree-sha1 = "91392e4664e82b3aa026086e35736047625f2a3f"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "6.52.0"
+version = "6.54.0"
 
     [deps.Symbolics.extensions]
     SymbolicsD3TreesExt = "D3Trees"
@@ -2595,9 +2619,9 @@ version = "0.2.0"
 
 [[deps.VectorizationBase]]
 deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
-git-tree-sha1 = "4ab62a49f1d8d9548a1c8d1a75e5f55cf196f64e"
+git-tree-sha1 = "d1d9a935a26c475ebffd54e9c7ad11627c43ea85"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.21.71"
+version = "0.21.72"
 
 [[deps.VersionParsing]]
 git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868"
diff --git a/benchmarks/NonStiffODE/Pleiades_wpd.jmd b/benchmarks/NonStiffODE/Pleiades_wpd.jmd
index f832f858e..97e8f7e75 100644
--- a/benchmarks/NonStiffODE/Pleiades_wpd.jmd
+++ b/benchmarks/NonStiffODE/Pleiades_wpd.jmd
@@ -4,7 +4,7 @@ author: Chris Rackauckas
 ---
 
 ```julia
-using OrdinaryDiffEq, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools
+using OrdinaryDiffEq, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools, IRKGaussLegendre
 
 f = (du, u, p, t) -> begin
     @inbounds begin
@@ -60,7 +60,9 @@ ODE.jl had to be discarded. The error estimate is off since it throws errors and
 setups = [Dict(:alg=>DP5())
           Dict(:alg=>dopri5())
           Dict(:alg=>Tsit5())
-          Dict(:alg=>Vern6())]
+          Dict(:alg=>Vern6())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol,
     save_everystep = false, numruns = 100, maxiters = 10000)
 plot(wp)
@@ -71,7 +73,9 @@ plot(wp)
 ```julia
 setups = [Dict(:alg=>DP5())
           Dict(:alg=>Tsit5())
-          Dict(:alg=>Vern6())]
+          Dict(:alg=>Vern6())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol, numruns = 100,
     maxiters = 10000, error_estimate = :L2, dense_errors = true)
 plot(wp)
@@ -93,6 +97,8 @@ setups = [Dict(:alg=>DP8())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>dop853())]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol,
     save_everystep = false, numruns = 100, maxiters = 1000)
@@ -106,6 +112,8 @@ setups = [Dict(:alg=>odex())
           Dict(:alg=>lsoda())
           Dict(:alg=>Vern6())
           Dict(:alg=>Tsit5())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>ddeabm())
           Dict(:alg=>ARKODE(Sundials.Explicit(), order = 6))]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol,
@@ -120,7 +128,9 @@ setups = [Dict(:alg=>DP8())
           Dict(:alg=>Vern6())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
-          Dict(:alg=>Vern9())]
+          Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol, numruns = 100,
     maxiters = 1000, error_estimate = :L2, dense_errors = true)
 plot(wp)
@@ -135,13 +145,15 @@ Adams-Bashforth-Moulton:
 setups = [Dict(:alg=>Tsit5())
           Dict(:alg=>Vern9())
           Dict(:alg=>VCABM())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>AitkenNeville(min_order = 1, max_order = 9, init_order = 4, threading = true))
           Dict(:alg=>ExtrapolationMidpointDeuflhard(
               min_order = 1, max_order = 9, init_order = 4, threading = true))
           Dict(:alg=>ExtrapolationMidpointHairerWanner(
               min_order = 2, max_order = 11, init_order = 4, threading = true))]
-solnames = ["Tsit5", "Vern9", "VCABM", "AitkenNeville",
-    "Midpoint Deuflhard", "Midpoint Hairer Wanner"]
+solnames = ["Tsit5", "Vern9", "VCABM", "IRKGL16", "IRKGL16 SIMD",
+    "AitkenNeville", "Midpoint Deuflhard", "Midpoint Hairer Wanner"]
 wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol = test_sol, names = solnames,
     save_everystep = false, verbose = false, numruns = 100)
 plot(wp)
diff --git a/benchmarks/NonStiffODE/Project.toml b/benchmarks/NonStiffODE/Project.toml
index 0732e4b5f..42e88a5fc 100644
--- a/benchmarks/NonStiffODE/Project.toml
+++ b/benchmarks/NonStiffODE/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+IRKGaussLegendre = "58bc7355-f626-4c51-96f2-1f8a038f95a2"
 LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
diff --git a/benchmarks/NonStiffODE/RigidBody_wpd.jmd b/benchmarks/NonStiffODE/RigidBody_wpd.jmd
index b9814b637..c2b6cfe87 100644
--- a/benchmarks/NonStiffODE/RigidBody_wpd.jmd
+++ b/benchmarks/NonStiffODE/RigidBody_wpd.jmd
@@ -5,7 +5,7 @@ author: Chris Rackauckas
 
 ```julia
 using OrdinaryDiffEq, ParameterizedFunctions, ODEInterfaceDiffEq, LSODA,
-      Sundials, DiffEqDevTools, StaticArrays, OrdinaryDiffEqSIMDRK
+      Sundials, DiffEqDevTools, StaticArrays, OrdinaryDiffEqSIMDRK, IRKGaussLegendre
 
 k(t) = 0.25*sin(t)^2
 
@@ -39,6 +39,8 @@ setups = [Dict(:alg=>DP5())
           Dict(:alg=>dopri5())
           Dict(:alg=>Tsit5())
           Dict(:alg=>Vern6())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Tsit5(), :prob_choice => 2)
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>MER5v2(), :prob_choice => 2)
@@ -61,6 +63,8 @@ setups = [Dict(:alg=>DP8())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>Vern7(), :prob_choice => 2)
           Dict(:alg=>Vern8(), :prob_choice => 2)
diff --git a/benchmarks/NonStiffODE/ThreeBody_wpd.jmd b/benchmarks/NonStiffODE/ThreeBody_wpd.jmd
index ee5839aa9..cd1672c55 100644
--- a/benchmarks/NonStiffODE/ThreeBody_wpd.jmd
+++ b/benchmarks/NonStiffODE/ThreeBody_wpd.jmd
@@ -4,7 +4,8 @@ author: Chris Rackauckas
 ---
 
 ```julia
-using OrdinaryDiffEq, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools, StaticArrays
+using OrdinaryDiffEq, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools, StaticArrays,
+      IRKGaussLegendre
 using Plots;
 gr()
 
@@ -80,6 +81,7 @@ setups = [Dict(:alg=>DP5())
           #Dict(:alg=>ode45()) #fails
           Dict(:alg=>BS5())
           Dict(:alg=>Tsit5())
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Tsit5(), :prob_choice => 2)
           Dict(:alg=>dopri5())];
 wp = WorkPrecisionSet(probs, abstols, reltols, setups; appxsol = test_sol,
@@ -94,7 +96,9 @@ setups = [Dict(:alg=>DP5(), :dense=>false)
           #Dict(:alg=>ode45()) # Fails
           Dict(:alg=>BS5(), :dense=>false)
           Dict(:alg=>Tsit5(), :dense=>false)
+          Dict(:alg=>RadauIIA5(), :dense=>false)
           Dict(:alg=>Tsit5(), :dense=>false, :prob_choice => 2)
+          Dict(:alg=>RadauIIA5(), :dense=>false, :prob_choice => 2)
           Dict(:alg=>dopri5())];
 wp = WorkPrecisionSet(probs, abstols, reltols, setups; appxsol = test_sol, numruns = 100)
 plot(wp)
@@ -107,6 +111,7 @@ setups = [Dict(:alg=>DP5())
           #Dict(:alg=>ode45()) #fails
           Dict(:alg=>BS5())
           Dict(:alg=>Tsit5())
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Tsit5(), :prob_choice => 2)
           Dict(:alg=>dopri5())];
 wp = WorkPrecisionSet(probs, abstols, reltols, setups; appxsol = test_sol, numruns = 100)
@@ -126,6 +131,8 @@ setups = [Dict(:alg=>DP5())
           Dict(:alg=>Vern7())
           Dict(:alg=>Vern8())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern6(), :prob_choice => 2)
           Dict(:alg=>Vern7(), :prob_choice => 2)
           Dict(:alg=>Vern8(), :prob_choice => 2)
@@ -163,6 +170,8 @@ Once again we separate ODE.jl because it fails. We also separate Sundials' `CVOD
 setups = [Dict(:alg=>DP5())
           Dict(:alg=>lsoda())
           Dict(:alg=>Vern8())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern8(), :prob_choice => 2)
           Dict(:alg=>ddeabm())
           Dict(:alg=>odex())
@@ -184,14 +193,17 @@ abstols = 1.0 ./ 10.0 .^ (3:13);
 reltols = 1.0 ./ 10.0 .^ (0:10);
 setups = [Dict(:alg=>Tsit5())
           Dict(:alg=>Vern9())
+          Dict(:alg=>IRKGL16(simd = false))
+          Dict(:alg=>IRKGL16(simd = true))
           Dict(:alg=>Vern9(), :prob_choice => 2)
           Dict(:alg=>AitkenNeville(min_order = 1, max_order = 9, init_order = 4, threading = true))
           Dict(:alg=>ExtrapolationMidpointDeuflhard(
               min_order = 1, max_order = 9, init_order = 4, threading = true))
           Dict(:alg=>ExtrapolationMidpointHairerWanner(
               min_order = 2, max_order = 11, init_order = 4, threading = true))]
-solnames = ["Tsit5", "Vern9", "Vern9 Static", "AitkenNeville",
-    "Midpoint Deuflhard", "Midpoint Hairer Wanner"]
+solnames = ["Tsit5", "Vern9", "IRKGL16", "IRKGL16 SIMD",
+    "Vern9 Static",
+    "AitkenNeville", "Midpoint Deuflhard", "Midpoint Hairer Wanner"]
 wp = WorkPrecisionSet(
     probs, abstols, reltols, setups; appxsol = test_sol, names = solnames,
     save_everystep = false, verbose = false, numruns = 100)