diff --git a/benchmarks/NonStiffBVP/Manifest.toml b/benchmarks/NonStiffBVP/Manifest.toml
index e53021083..072246d8f 100644
--- a/benchmarks/NonStiffBVP/Manifest.toml
+++ b/benchmarks/NonStiffBVP/Manifest.toml
@@ -2111,4 +2111,4 @@ version = "3.5.0+0"
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
 git-tree-sha1 = "9c304562909ab2bab0262639bd4f444d7bc2be37"
 uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
-version = "1.4.1+1"
+version = "1.4.1+1"
\ No newline at end of file
diff --git a/benchmarks/NonStiffBVP/Project.toml b/benchmarks/NonStiffBVP/Project.toml
index afdc20a7b..d3c76fdf9 100644
--- a/benchmarks/NonStiffBVP/Project.toml
+++ b/benchmarks/NonStiffBVP/Project.toml
@@ -16,4 +16,4 @@ DiffEqDevTools = "2.42.0"
 ODEInterface = "0.5"
 OrdinaryDiffEq = "6"
 Plots = "1.39"
-SciMLBenchmarks = "0.1"
+SciMLBenchmarks = "0.1"
\ No newline at end of file
diff --git a/benchmarks/NonStiffBVP/nonlinear_wpi.jmd b/benchmarks/NonStiffBVP/nonlinear_wpi.jmd
new file mode 100644
index 000000000..befb0bc04
--- /dev/null
+++ b/benchmarks/NonStiffBVP/nonlinear_wpi.jmd
@@ -0,0 +1,247 @@
+
+This is a benchmark of MIRK methods and Shooting method on some nonlinear boundary value problems. The testing BVPs are a set of standard BVP test problems as described [here](https://archimede.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=29)
+
+# Setup
+
+Fetch required packages.
+
+```julia
+using BoundaryValueDiffEq, OrdinaryDiffEq, DiffEqDevTools, BenchmarkTools, BVProblemLibrary, Plots
+```
+
+Set up the benchmarked solvers.
+
+```julia
+setups = [ Dict(:alg=>MIRK2(), :dt=>0.1),
+            Dict(:alg=>MIRK3(), :dt=>0.1),
+            Dict(:alg=>MIRK4(), :dt=>0.1),
+            Dict(:alg=>MIRK5(), :dt=>0.1),
+            Dict(:alg=>MIRK6(), :dt=>0.1),
+            Dict(:alg=>Shooting(Tsit5())),
+            Dict(:alg=>Shooting(Vern7())),
+            Dict(:alg=>MultipleShooting(10,Tsit5())),
+            Dict(:alg=>MultipleShooting(10,Vern7()))]
+labels = ["MIRK2";
+               "MIRK3";
+               "MIRK4";
+               "MIRK5";
+               "MIRK6";
+               "Shooting(Tsit5)";
+               "Shooting(Vern7)";
+               "MultipleShooting(Tsit5)";
+               "MultipleShooting(Vern7)"];
+```
+
+Sets tolerances.
+
+```julia
+abstols = 1.0 ./ 10.0 .^ (1:4)
+reltols = 1.0 ./ 10.0 .^ (1:4);
+```
+
+Function that checks retcode.
+
+```julia
+function retcode_check(setups, labels, prob)
+    success_setups = []
+    success_labels = []
+    for (setup, label) in zip(setups, labels)
+        temp_sol = solve(prob, setup[:alg], abstol=1/10^8,reltol=1/10^8, dt=0.1)
+        if SciMLBase.successful_retcode(temp_sol.retcode)
+            push!(success_setups, setup)
+            push!(success_labels, label)
+        end
+    end
+    return success_setups, success_labels
+end
+```
+
+# Benchmarks
+
+We here run the BVP benchmarks for all the MIRK methods.
+
+## Nonlinear boundary value problems
+
+### Nonlinear BVP 1
+
+```julia
+prob_1 = BVProblemLibrary.prob_bvp_nonlinear_1
+sol_1 = solve(prob_1, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_1 = TestSolution(sol_1.t, sol_1.u)
+setups_1, labels_1 = retcode_check(setups, labels, prob_1)
+wp_1 = WorkPrecisionSet(prob_1, abstols, reltols, setups_1; names = labels_1, print_names = true, appxsol = testsol_1, 
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_1, title="Nonlinear BVP 1")
+```
+
+### Nonlinear BVP 2
+
+```julia
+prob_2 = BVProblemLibrary.prob_bvp_nonlinear_2
+sol_2 = solve(prob_2, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_2 = TestSolution(sol_2.t, sol_2.u)
+setups_2, labels_2 = retcode_check(setups, labels, prob_2)
+wp_2 = WorkPrecisionSet(prob_2, abstols, reltols, setups_2; names = labels_2, print_names = true, appxsol = testsol_2,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_2, title="Nonlinear BVP 2")
+```
+
+### Nonlinear BVP 3
+
+```julia
+prob_3 = BVProblemLibrary.prob_bvp_nonlinear_3
+sol_3 = solve(prob_3, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_3 = TestSolution(sol_3.t, sol_3.u)
+setups_3, labels_3 = retcode_check(setups, labels, prob_3)
+wp_3 = WorkPrecisionSet(prob_3, abstols, reltols, setups_3; names = labels_3, print_names = true, appxsol = testsol_3,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_3, title="Nonlinear BVP 3")
+```
+
+### Noninear BVP 4
+
+```julia
+prob_4 = BVProblemLibrary.prob_bvp_nonlinear_4
+sol_4 = solve(prob_4, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_4 = TestSolution(sol_4.t, sol_4.u)
+setups_4, labels_4 = retcode_check(setups, labels, prob_4)
+wp_4 = WorkPrecisionSet(prob_4, abstols, reltols, setups_4; names = labels_4, print_names = true, appxsol = testsol_4,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_4, title="Nonlinear BVP 4")
+```
+
+### Nonlinear BVP 5
+
+```julia
+prob_5 = BVProblemLibrary.prob_bvp_nonlinear_5
+sol_5 = solve(prob_5, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_5 = TestSolution(sol_5.t, sol_5.u)
+setups_5, labels_5 = retcode_check(setups, labels, prob_5)
+wp_5 = WorkPrecisionSet(prob_5, abstols, reltols, setups_5; names = labels_5, print_names = true, appxsol = testsol_5,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_5, title="Nonlinear BVP 5")
+```
+
+### Nonlinear BVP 6
+
+### ```julia
+### prob_6 = BVProblemLibrary.prob_bvp_nonlinear_6
+### sol_6 = solve(prob_6, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+### testsol_6 = TestSolution(sol_6.t, sol_6.u)
+### wp_6 = WorkPrecisionSet(prob_6, abstols, reltols, setups_6; names = labels_6, print_names = true, appxsol = testsol_6,
+###                         dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+### plot(wp_6)
+### ```
+
+### Nonlinear BVP 7
+
+```julia
+prob_7 = BVProblemLibrary.prob_bvp_nonlinear_7
+sol_7 = solve(prob_7, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_7 = TestSolution(sol_7.t, sol_7.u)
+setups_7, labels_7 = retcode_check(setups, labels, prob_7)
+wp_7 = WorkPrecisionSet(prob_7, abstols, reltols, setups_7; names = labels_7, print_names = true, appxsol = testsol_7,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_7, title="Nonlinear BVP 7")
+```
+
+### Nonlinear BVP 8
+
+```julia
+prob_8 = BVProblemLibrary.prob_bvp_nonlinear_8
+sol_8 = solve(prob_8, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_8 = TestSolution(sol_8.t, sol_8.u)
+setups_8, labels_8 = retcode_check(setups, labels, prob_8)
+wp_8 = WorkPrecisionSet(prob_8, abstols, reltols, setups_8; names = labels_8, print_names = true, appxsol = testsol_8,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_8, title="Nonlinear BVP 8")
+```
+
+### Nonlinear BVP 9
+
+```julia
+prob_9 = BVProblemLibrary.prob_bvp_nonlinear_9
+sol_9 = solve(prob_9, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_9 = TestSolution(sol_9.t, sol_9.u)
+setups_9, labels_9 = retcode_check(setups, labels, prob_9)
+wp_9 = WorkPrecisionSet(prob_9, abstols, reltols, setups_9; names = labels_9, print_names = true, appxsol = testsol_9,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_9, title="Nonlinear BVP 9")
+```
+
+### Nonlinear BVP 10
+
+```julia
+prob_10 = BVProblemLibrary.prob_bvp_nonlinear_10
+sol_10 = solve(prob_10, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_10 = TestSolution(sol_10.t, sol_10.u)
+setups_10, labels_10 = retcode_check(setups, labels, prob_10)
+wp_10 = WorkPrecisionSet(prob_10, abstols, reltols, setups_10; names = labels_10, print_names = true, appxsol = testsol_10,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_10, title="Nonlinear BVP 10")
+```
+
+### Nonlinear BVP 11
+
+```julia
+prob_11 = BVProblemLibrary.prob_bvp_nonlinear_11
+sol_11 = solve(prob_11, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_11 = TestSolution(sol_11.t, sol_11.u)
+setups_11, labels_11 = retcode_check(setups, labels, prob_11)
+wp_11 = WorkPrecisionSet(prob_11, abstols, reltols, setups_11; names = labels_11, print_names = true, appxsol = testsol_11,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_11, title="Nonlinear BVP 11")
+```
+
+### Nonlinear BVP 12
+
+```julia
+prob_12 = BVProblemLibrary.prob_bvp_nonlinear_12
+sol_12 = solve(prob_12, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_12 = TestSolution(sol_12.t, sol_12.u)
+setups_12, labels_12 = retcode_check(setups, labels, prob_12)
+wp_12 = WorkPrecisionSet(prob_12, abstols, reltols, setups_12; names = labels_12, print_names = true, appxsol = testsol_12,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_12, title="Nonlinear BVP 12")
+```
+
+### Nonlinear BVP 13
+
+```julia
+prob_13 = BVProblemLibrary.prob_bvp_nonlinear_13
+sol_13 = solve(prob_13, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_13 = TestSolution(sol_13.t, sol_13.u)
+setups_13, labels_13 = retcode_check(setups, labels, prob_13)
+wp_13 = WorkPrecisionSet(prob_13, abstols, reltols, setups_13; names = labels_13, print_names = true, appxsol = testsol_13,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_13, title="Nonlinear BVP 13")
+```
+
+### Nonlinear BVP 14
+
+```julia
+prob_14 = BVProblemLibrary.prob_bvp_nonlinear_14
+sol_14 = solve(prob_14, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_14 = TestSolution(sol_14.t, sol_14.u)
+setups_14, labels_14 = retcode_check(setups, labels, prob_14)
+wp_14 = WorkPrecisionSet(prob_14, abstols, reltols, setups_14; names = labels_14, print_names = true, appxsol = testsol_14,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000, timeseries_errors=false,verbose=false)
+plot(wp_14, title="Nonlinear BVP 14")
+```
+
+### Nonlinear BVP 15
+
+```julia
+prob_15 = BVProblemLibrary.prob_bvp_nonlinear_15
+sol_15 = solve(prob_15, MIRK3(), abstol=1e-8, reltol=1e-8, dt = 0.1)
+testsol_15 = TestSolution(sol_15.t, sol_15.u)
+setups_15, labels_15 = retcode_check(setups, labels, prob_15)
+wp_15 = WorkPrecisionSet(prob_15, abstols, reltols, setups_15; names = labels_15, print_names = true, appxsol = testsol_15,
+                        dense=false, save_everystep=false,numruns=10,maxiters=1000,timeseries_errors=false,verbose=false)
+plot(wp_15, title="Nonlinear BVP 15")
+```
+
+# ```julia, echo = false
+# using SciMLBenchmarks
+# SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
+# ```
\ No newline at end of file