diff --git a/lib/OrdinaryDiffEqLowStorageRK/Project.toml b/lib/OrdinaryDiffEqLowStorageRK/Project.toml
index 5e86e1345e..4c0497a761 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/Project.toml
+++ b/lib/OrdinaryDiffEqLowStorageRK/Project.toml
@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqLowStorageRK"
 uuid = "b0944070-b475-4768-8dec-fb6eb410534d"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "3.3.0"
+version = "3.4.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
index 0495abd928..e437c6c346 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
@@ -60,6 +60,7 @@ alg_order(alg::SHLDDRK_2N) = 4
 alg_order(alg::SHLDDRK52) = 2
 
 isfsal(alg::ORK256) = false
+isfsal(alg::RK46NL) = false
 isfsal(alg::CarpenterKennedy2N54) = false
 isfsal(alg::DGLDDRK84_F) = false
 isfsal(alg::DGLDDRK73_C) = false
@@ -73,6 +74,7 @@ isfsal(alg::NDBLSRK124) = false
 isfsal(alg::NDBLSRK144) = false
 
 uses_uprev(alg::ORK256, adaptive::Bool) = false
+uses_uprev(alg::RK46NL, adaptive::Bool) = false
 uses_uprev(alg::SHLDDRK64, adaptive::Bool) = false
 uses_uprev(alg::CarpenterKennedy2N54, adaptive::Bool) = false
 uses_uprev(alg::NDBLSRK124, adaptive::Bool) = false
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl
index 48a84d8274..a2946c6635 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl
@@ -195,12 +195,16 @@ end
 @doc explicit_rk_docstring(
     "6-stage, fourth order low-stage, low-dissipation, low-dispersion scheme.
 Fixed timestep only.", "RK46NL",
-    references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003"
+    references = "Julien Berland, Christophe Bogey, Christophe Bailly. Low-Dissipation and Low-Dispersion Fourth-Order Runge-Kutta Algorithm. Computers & Fluids, 35(10), pp 1459-1463, 2006. doi: https://doi.org/10.1016/j.compfluid.2005.04.003",
+    extra_keyword_description = """- `williamson_condition`: allows for an optimization that allows fusing broadcast expressions with the function call `f`. However, it only works for `Array` types.
+    """,
+    extra_keyword_default = "williamson_condition = false"
 )
 Base.@kwdef struct RK46NL{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm
     stage_limiter!::StageLimiter = trivial_limiter!
     step_limiter!::StepLimiter = trivial_limiter!
     thread::Thread = Serial()
+    williamson_condition::Bool = false
 end
 
 @doc explicit_rk_docstring(
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl
index 6a436b6ce6..0c71fb3aee 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl
@@ -47,81 +47,30 @@ function ORK256ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
     return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
 end
 
-function alg_cache(
-        alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{false}, verbose
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    return RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
-end
+function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
+    A2 = convert(T, -0.737101392796)
+    A3 = convert(T, -1.634740794343)
+    A4 = convert(T, -0.74473900378)
+    A5 = convert(T, -1.469897351522)
+    A6 = convert(T, -2.813971388035)
+    A2end = (A2, A3, A4, A5, A6)
 
-struct RK46NLConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
-    α2::T
-    α3::T
-    α4::T
-    α5::T
-    α6::T
-    β1::T
-    β2::T
-    β3::T
-    β4::T
-    β5::T
-    β6::T
-    c2::T2
-    c3::T2
-    c4::T2
-    c5::T2
-    c6::T2
-
-    function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
-        α2 = T(-0.737101392796)
-        α3 = T(-1.634740794343)
-        α4 = T(-0.74473900378)
-        α5 = T(-1.469897351522)
-        α6 = T(-2.813971388035)
-        β1 = T(0.032918605146)
-        β2 = T(0.8232569982)
-        β3 = T(0.3815309489)
-        β4 = T(0.200092213184)
-        β5 = T(1.718581042715)
-        β6 = T(0.27)
-        c2 = T2(0.032918605146)
-        c3 = T2(0.249351723343)
-        c4 = T2(0.466911705055)
-        c5 = T2(0.582030414044)
-        c6 = T2(0.847252983783)
-        return new{T, T2}(α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6)
-    end
-end
+    B1 = convert(T, 0.032918605146)
+    B2 = convert(T, 0.8232569982)
+    B3 = convert(T, 0.3815309489)
+    B4 = convert(T, 0.200092213184)
+    B5 = convert(T, 1.718581042715)
+    B6 = convert(T, 0.27)
+    B2end = (B2, B3, B4, B5, B6)
 
-function alg_cache(
-        alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits},
-        ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
-        dt, reltol, p, calck,
-        ::Val{true}, verbose
-    ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    tmp = zero(u)
-    k = zero(rate_prototype)
-    fsalfirst = zero(rate_prototype)
-    tab = RK46NLConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits))
-    return RK46NLCache(
-        u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!,
-        alg.thread
-    )
-end
+    c2 = convert(T2, 0.032918605146)
+    c3 = convert(T2, 0.249351723343)
+    c4 = convert(T2, 0.466911705055)
+    c5 = convert(T2, 0.582030414044)
+    c6 = convert(T2, 0.847252983783)
+    c2end = (c2, c3, c4, c5, c6)
 
-@cache struct RK46NLCache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <:
-    LowStorageRKMutableCache
-    u::uType
-    uprev::uType
-    k::rateType
-    tmp::uType
-    fsalfirst::rateType
-    tab::TabType
-    stage_limiter!::StageLimiter
-    step_limiter!::StepLimiter
-    thread::Thread
+    return LowStorageRK2NConstantCache(A2end, B1, B2end, c2end)
 end
 
 function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T, T2}
@@ -4357,11 +4306,12 @@ end
 
 # --- LowStorageRK2N ---
 const LowStorageRK2NAlgorithm = Union{
-    ORK256, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F,
+    ORK256, RK46NL, CarpenterKennedy2N54, SHLDDRK64, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F,
     NDBLSRK124, NDBLSRK134, NDBLSRK144,
 }
 
 _lowstorage_2n_tableau(::ORK256, ::Type{T}, ::Type{T2}) where {T, T2} = ORK256ConstantCache(T, T2)
+_lowstorage_2n_tableau(::RK46NL, ::Type{T}, ::Type{T2}) where {T, T2} = RK46NLConstantCache(T, T2)
 _lowstorage_2n_tableau(::CarpenterKennedy2N54, ::Type{T}, ::Type{T2}) where {T, T2} = CarpenterKennedy2N54ConstantCache(T, T2)
 _lowstorage_2n_tableau(::SHLDDRK64, ::Type{T}, ::Type{T2}) where {T, T2} = SHLDDRK64ConstantCache(T, T2)
 _lowstorage_2n_tableau(::DGLDDRK73_C, ::Type{T}, ::Type{T2}) where {T, T2} = DGLDDRK73_CConstantCache(T, T2)
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl
index 24ef9f39fa..21cf4e024b 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl
@@ -903,96 +903,6 @@ end
     OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
 end
 
-function initialize!(integrator, cache::RK46NLCache)
-    (; k, fsalfirst) = cache
-
-    integrator.kshortsize = 1
-    resize!(integrator.k, integrator.kshortsize)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation
-    return OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-end
-
-@muladd function perform_step!(integrator, cache::RK46NLCache, repeat_step = false)
-    (; t, dt, uprev, u, f, p) = integrator
-    (; k, fsalfirst, tmp, thread) = cache
-    stage_limiter! = integrator.opts.stage_limiter!
-    (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache.tab
-
-    # u1
-    @.. broadcast = false thread = thread tmp = dt * fsalfirst
-    @.. broadcast = false thread = thread u = uprev + β1 * tmp
-    stage_limiter!(u, integrator, p, t + c2 * dt)
-    # u2
-    f(k, u, p, t + c2 * dt)
-    @.. broadcast = false thread = thread tmp = α2 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β2 * tmp
-    stage_limiter!(u, integrator, p, t + c3 * dt)
-    # u3
-    f(k, u, p, t + c3 * dt)
-    @.. broadcast = false thread = thread tmp = α3 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β3 * tmp
-    stage_limiter!(u, integrator, p, t + c4 * dt)
-    # u4
-    f(k, u, p, t + c4 * dt)
-    @.. broadcast = false thread = thread tmp = α4 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β4 * tmp
-    stage_limiter!(u, integrator, p, t + c5 * dt)
-    # u5 = u
-    f(k, u, p, t + c5 * dt)
-    @.. broadcast = false thread = thread tmp = α5 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β5 * tmp
-    stage_limiter!(u, integrator, p, t + c6 * dt)
-
-    f(k, u, p, t + c6 * dt)
-    @.. broadcast = false thread = thread tmp = α6 * tmp + dt * k
-    @.. broadcast = false thread = thread u = u + β6 * tmp
-    stage_limiter!(u, integrator, p, t + dt)
-
-    f(k, u, p, t + dt)
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6)
-end
-
-function initialize!(integrator, cache::RK46NLConstantCache)
-    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
-    integrator.kshortsize = 1
-    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
-
-    # Avoid undefined entries if k is an array of arrays
-    integrator.fsallast = zero(integrator.fsalfirst)
-    return integrator.k[1] = integrator.fsalfirst
-end
-
-@muladd function perform_step!(integrator, cache::RK46NLConstantCache, repeat_step = false)
-    (; t, dt, uprev, u, f, p) = integrator
-    (; α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6) = cache
-
-    # u1
-    tmp = dt * integrator.fsalfirst
-    u = uprev + β1 * tmp
-    # u2
-    tmp = α2 * tmp + dt * f(u, p, t + c2 * dt)
-    u = u + β2 * tmp
-    # u3
-    tmp = α3 * tmp + dt * f(u, p, t + c3 * dt)
-    u = u + β3 * tmp
-    # u4
-    tmp = α4 * tmp + dt * f(u, p, t + c4 * dt)
-    u = u + β4 * tmp
-    # u5 = u
-    tmp = α5 * tmp + dt * f(u, p, t + c5 * dt)
-    u = u + β5 * tmp
-    # u6
-    tmp = α6 * tmp + dt * f(u, p, t + c6 * dt)
-    u = u + β6 * tmp
-
-    integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd
-    OrdinaryDiffEqCore.increment_nf!(integrator.stats, 6)
-    integrator.k[1] = integrator.fsalfirst
-    integrator.u = u
-end
-
 function initialize!(integrator, cache::SHLDDRK52ConstantCache)
     integrator.kshortsize = 2
     integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
diff --git a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl
index 27724b5889..df998c58da 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl
@@ -143,6 +143,64 @@ end
     @test sol_old.u[end] ≈ sol_new.u[end]
 end
 
+@testset "RK46NL" begin
+    alg = RK46NL(; williamson_condition = true)
+    alg2 = RK46NL(; williamson_condition = false)
+    dts = 1 ./ 2 .^ (7:-1:3)
+    for prob in test_problems_only_time
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    for prob in test_problems_linear
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    for prob in test_problems_nonlinear
+        sim = test_convergence(dts, prob, alg)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+        sim = test_convergence(dts, prob, alg2)
+        @test sim.𝒪est[:final] ≈ OrdinaryDiffEqLowStorageRK.alg_order(alg) atol = testTol
+    end
+    integ = init(
+        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
+    integ = init(
+        prob_ode_large, alg, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 2
+    integ = init(
+        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 4
+    integ = init(
+        prob_ode_large, alg2, dt = 1.0e-2, save_start = false, save_end = false,
+        save_everystep = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 3
+    # test whether aliasing u0 is bad
+    new_prob_ode_nonlinear_inplace = ODEProblem(
+        prob_ode_nonlinear_inplace.f, [1.0],
+        (0.0, 0.5)
+    )
+    sol_old = solve(
+        prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
+        save_start = false
+    )
+    sol_new = solve(
+        new_prob_ode_nonlinear_inplace, alg, dt = 1.0e-4, save_everystep = false,
+        save_start = false, alias = ODEAliasSpecifier(alias_u0 = true)
+    )
+    @test sol_old.u[end] ≈ sol_new.u[end]
+end
+
 @testset "CarpenterKennedy2N54" begin
     alg = CarpenterKennedy2N54(; williamson_condition = true)
     alg2 = CarpenterKennedy2N54(; williamson_condition = false)