diff --git a/docs/Project.toml b/docs/Project.toml
index 573533e21..d87833c7d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -9,6 +9,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
diff --git a/docs/pages.jl b/docs/pages.jl
index e9999c82c..8552f3d19 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -6,7 +6,8 @@ pages = ["index.md",
         "tutorials/large_systems.md",
         "tutorials/modelingtoolkit.md",
         "tutorials/small_compile.md",
-        "tutorials/iterator_interface.md"],
+        "tutorials/iterator_interface.md",
+        "tutorials/optimizing_parameterized_ode.md",],
     "Basics" => Any["basics/NonlinearProblem.md",
         "basics/NonlinearFunctions.md",
         "basics/solve.md",
diff --git a/docs/src/tutorials/optimizing_parameterized_ode.md b/docs/src/tutorials/optimizing_parameterized_ode.md
new file mode 100644
index 000000000..44ba65eb5
--- /dev/null
+++ b/docs/src/tutorials/optimizing_parameterized_ode.md
@@ -0,0 +1,76 @@
+# [Optimizing a Parameterized ODE](@id optimizing-parameterized-ode)
+
+Let us fit a parameterized ODE to some data. We will use the Lotka-Volterra model as an
+example. We will use Single Shooting to fit the parameters.
+
+```@example parameterized_ode
+using OrdinaryDiffEq, NonlinearSolve, Plots
+```
+
+Let us simulate some real data from the Lotka-Volterra model.
+
+```@example parameterized_ode
+function lotka_volterra!(du, u, p, t)
+    x, y = u
+    α, β, δ, γ = p
+    du[1] = dx = α * x - β * x * y
+    du[2] = dy = -δ * y + γ * x * y
+end
+
+# Initial condition
+u0 = [1.0, 1.0]
+
+# Simulation interval and intermediary points
+tspan = (0.0, 2.0)
+tsteps = 0.0:0.1:10.0
+
+# LV equation parameter. p = [α, β, δ, γ]
+p = [1.5, 1.0, 3.0, 1.0]
+
+# Setup the ODE problem, then solve
+prob = ODEProblem(lotka_volterra!, u0, tspan, p)
+sol = solve(prob, Tsit5(); saveat = tsteps)
+
+# Plot the solution
+using Plots
+plot(sol; linewidth = 3)
+savefig("LV_ode.png")
+```
+
+Let us now formulate the parameter estimation as a Nonlinear Least Squares Problem.
+
+```@example parameterized_ode
+function loss_function(ode_param, data)
+    sol = solve(prob, Tsit5(); p = ode_param, saveat = tsteps)
+    return vec(reduce(hcat, sol.u)) .- data
+end
+
+p_init = zeros(4)
+
+nlls_prob = NonlinearLeastSquaresProblem(loss_function, p_init, vec(reduce(hcat, sol.u)))
+```
+
+Now, we can use any NLLS solver to solve this problem.
+
+```@example parameterized_ode
+res = solve(nlls_prob, LevenbergMarquardt(); maxiters = 1000, show_trace = Val(true),
+    trace_level = TraceAll())
+```
+
+We can also use Trust Region methods.
+
+```@example parameterized_ode
+res = solve(nlls_prob, TrustRegion(); maxiters = 1000, show_trace = Val(true),
+    trace_level = TraceAll())
+```
+
+Let's plot the solution.
+
+```@example parameterized_ode
+prob2 = remake(prob; tspan = (0.0, 10.0))
+sol_fit = solve(prob2, Tsit5(); p = res.u)
+sol_true = solve(prob2, Tsit5(); p = p)
+plot(sol_true; linewidth = 3)
+plot!(sol_fit; linewidth = 3, linestyle = :dash)
+savefig("LV_ode_fit.png")
+```