diff --git a/test/Project.toml b/test/Project.toml
index c022d8b8e..7dd42e1dc 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,10 +3,10 @@ CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 MadNLPMumps = "3b83494e-c0a4-4895-918b-9157a7a085a1"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -16,10 +16,10 @@ CTModels = "0.6"
 DifferentiationInterface = "0.7"
 ForwardDiff = "0.10, 1.0"
 LinearAlgebra = "1"
-MINPACK = "1"
 MadNLP = "0.8"
 MadNLPMumps = "0.5"
 NLPModelsIpopt = "0.10"
+NonlinearSolve = "4"
 OrdinaryDiffEq = "6"
 SplitApplyCombine = "1"
 Test = "1"
diff --git a/test/indirect/test_goddard_indirect.jl b/test/indirect/test_goddard_indirect.jl
index efa5facf5..32c22a6b9 100644
--- a/test/indirect/test_goddard_indirect.jl
+++ b/test/indirect/test_goddard_indirect.jl
@@ -42,39 +42,20 @@ function test_goddard_indirect()
     end
 
     # tests
-    t1 = 0.025246759388000528
-    t2 = 0.061602092906721286
-    t3 = 0.10401664867856217
-    tf = 0.20298394547952422
-    p0 = [3.9428857983400074, 0.14628855388160236, 0.05412448008321635]
+    p0 = [3.9457646586891744, 0.15039559623165552, 0.05371271293970545]
+    t1 = 0.023509684041879215
+    t2 = 0.059737380899876
+    t3 = 0.10157134842432228
+    tf = 0.20204744057100849
 
-    # test shooting function
+    # test shooting function with solve from NonlinearSolve
     s = zeros(eltype(p0), 7)
-    shoot!(s, p0, t1, t2, t3, tf)
-    s_guess_sol = [
-        -0.02456074767656735,
-        -0.05699760226157302,
-        0.0018629693253921868,
-        -0.027013078908634858,
-        -0.21558816838342798,
-        -0.0121146739026253,
-        0.015713236406057297,
-    ]
-    @test s ≈ s_guess_sol atol = 1e-6
-
-    # solve and compare
     ξ0 = [p0; t1; t2; t3; tf]
-    backend = AutoForwardDiff()
-    nle! = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
-    jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ)
-    indirect_sol = fsolve(nle!, jnle!, ξ0; show_trace=true)
-
-    p0 = indirect_sol.x[1:3]
-    t1 = indirect_sol.x[4]
-    t2 = indirect_sol.x[5]
-    t3 = indirect_sol.x[6]
-    tf = indirect_sol.x[7]
-
-    shoot!(s, p0, t1, t2, t3, tf)
+    shoot!(s, ξ, λ) = shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
+    prob = NonlinearProblem(shoot!, ξ0)
+    sol = solve(prob)
+    ξ = sol.u
+    p0, t1, t2, t3, tf = ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7]
     @test norm(s) < 1e-6
+
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index a93e4a294..a2c007763 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,10 +5,10 @@ using MadNLP
 using MadNLPMumps
 using LinearAlgebra
 using OrdinaryDiffEq
-using MINPACK
 using DifferentiationInterface
 using ForwardDiff: ForwardDiff
 using SplitApplyCombine # for flatten in some tests
+using NonlinearSolve
 
 # NB some direct tests use functional definition and are `using CTModels`