diff --git a/lib/SimpleNonlinearSolve/Project.toml b/lib/SimpleNonlinearSolve/Project.toml
index f1608e1ce..53749a72f 100644
--- a/lib/SimpleNonlinearSolve/Project.toml
+++ b/lib/SimpleNonlinearSolve/Project.toml
@@ -1,7 +1,7 @@
 name = "SimpleNonlinearSolve"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 authors = ["SciML"]
-version = "1.0.2"
+version = "1.0.3"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/lib/SimpleNonlinearSolve/src/nlsolve/lbroyden.jl b/lib/SimpleNonlinearSolve/src/nlsolve/lbroyden.jl
index 6dcac6a37..a78611a0a 100644
--- a/lib/SimpleNonlinearSolve/src/nlsolve/lbroyden.jl
+++ b/lib/SimpleNonlinearSolve/src/nlsolve/lbroyden.jl
@@ -68,8 +68,8 @@ end
         d = dot(vᵀ, δf)
         @bb @. δx = (δx - mvec) / d
 
-        selectdim(U, 2, mod1(i, η)) .= δx
-        selectdim(Vᵀ, 1, mod1(i, η)) .= vᵀ
+        selectdim(U, 2, mod1(i, η)) .= _vec(δx)
+        selectdim(Vᵀ, 1, mod1(i, η)) .= _vec(vᵀ)
 
         _U = selectdim(U, 2, 1:min(η, i))
         _Vᵀ = selectdim(Vᵀ, 1, 1:min(η, i))
@@ -93,7 +93,7 @@ function _rmatvec!!(y, xᵀU, U, Vᵀ, x)
     x_ = vec(x)
     xᵀU_ = view(xᵀU, 1:η)
     @bb xᵀU_ = transpose(U) × x_
-    @bb y = transpose(Vᵀ) × xᵀU_
+    @bb y = transpose(Vᵀ) × vec(xᵀU_)
     @bb @. y -= x
     return y
 end
@@ -108,7 +108,7 @@ function _matvec!!(y, Vᵀx, U, Vᵀ, x)
     x_ = vec(x)
     Vᵀx_ = view(Vᵀx, 1:η)
     @bb Vᵀx_ = Vᵀ × x_
-    @bb y = U × Vᵀx_
+    @bb y = U × vec(Vᵀx_)
     @bb @. y -= x
     return y
 end
diff --git a/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl b/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
index 66f6a3d0c..455aac91a 100644
--- a/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
+++ b/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
@@ -1,12 +1,12 @@
 using SimpleNonlinearSolve
 
 ff(u, p) = u .* u .- p
-u0 = rand(2, 2)
+u0 = ones(2, 3)
 p = 2.0
 vecprob = NonlinearProblem(ff, vec(u0), p)
 prob = NonlinearProblem(ff, u0, p)
 
 @testset "$(alg)" for alg in (SimpleKlement(), SimpleBroyden(), SimpleNewtonRaphson(),
-    SimpleDFSane(), SimpleLimitedMemoryBroyden(), SimpleTrustRegion())
+    SimpleDFSane(), SimpleLimitedMemoryBroyden(; threshold = Val(2)), SimpleTrustRegion())
     @test vec(solve(prob, alg).u) ≈ solve(vecprob, alg).u
 end