diff --git a/Manifest.toml b/Manifest.toml
index 7dc10e82..f41e7b01 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.0"
 manifest_format = "2.0"
-project_hash = "4a03429a978da71126f1de364ece06ed1b1e7c3b"
+project_hash = "7e5807e059f83e7dc9033cceb2df5b3a68479e06"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245"
@@ -570,9 +570,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"]
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
-git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
+git-tree-sha1 = "73d1214fec245096717847c62d389a5d2ac86504"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.21.1"
+version = "2.22.0"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
diff --git a/Project.toml b/Project.toml
index 01855d50..f4bd63c0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,6 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/docs/figs/ch2-04.jl b/docs/figs/ch2-04.jl
index 259fb1b5..7643d4e7 100644
--- a/docs/figs/ch2-04.jl
+++ b/docs/figs/ch2-04.jl
@@ -15,4 +15,6 @@ for k in 1:3
     plot!(fig, t -> 3 * exp(-k*t), 0., 5., label = "exp(-$(k)t)")
 end
 
-plot!(fig, xlim = (0, 5), ylim=(0, 3.2), xlabel="Time", ylabel="Concentration") |> PNG
+plot!(fig, xlim = (0, 5), ylim=(0, 3.2), xlabel="Time", ylabel="Concentration")
+
+fig |> PNG
diff --git a/docs/figs/ch2-p246.jl b/docs/figs/ch2-p246.jl
index b1b47df8..944e5bd0 100644
--- a/docs/figs/ch2-p246.jl
+++ b/docs/figs/ch2-p246.jl
@@ -8,6 +8,4 @@ Plots.default(linewidth=2)
 PNG(fig) = display("image/png", fig)
 
 #---
-fig = ODEProblem((u, p, t) -> p * (1. - u), 0., 10., 1.) |> solve |> plot
-
-fig |> PNG
+ODEProblem((u, p, t) -> p * (1. - u), 0., 10., 1.) |> solve |> plot |> PNG
diff --git a/docs/figs/ch4-01.jl b/docs/figs/ch4-01.jl
index f87a2e79..06146eae 100644
--- a/docs/figs/ch4-01.jl
+++ b/docs/figs/ch4-01.jl
@@ -5,7 +5,6 @@ Steady states and phase plots in an asymmetric network.
 ===#
 
 using DifferentialEquations
-using LabelledArrays
 using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
@@ -41,11 +40,11 @@ end
 
 ps1 = (k1=20., k2=5., k3=5., k4=5., k5=2., n=4.)
 u0s = [
-    LVector(a=0.0, b=0.0),
-    LVector(a=0.5, b=0.6),
-    LVector(a=0.17, b=1.1),
-    LVector(a=0.25, b=1.9),
-    LVector(a=1.85, b=1.7),
+    [0.0, 0.0],
+    [0.5, 0.6],
+    [0.17, 1.1],
+    [0.25, 1.9],
+    [1.85, 1.7],
 ]
 
 tend = 1.5
diff --git a/docs/figs/ch4-07.jl b/docs/figs/ch4-07.jl
index fa9c7441..c74c0e5d 100644
--- a/docs/figs/ch4-07.jl
+++ b/docs/figs/ch4-07.jl
@@ -2,8 +2,6 @@
 # Symmetric (bistable) biological networks.
 
 using DifferentialEquations
-using LabelledArrays
-using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
@@ -17,13 +15,13 @@ hil(x, k, n) = hil(x^n, k^n)
 #---
 function ∂A47(u, p, t)
     a, b = u
-    @unpack k1, k2, k3, k4, n1, n2 = p
+    k1, k2, k3, k4, n1, n2 = p
     return da = k1 * hil(1, b, n1) - k3 * a
 end
 
 function ∂B47(u, p, t)
     a, b = u
-    @unpack k1, k2, k3, k4, n1, n2 = p
+    k1, k2, k3, k4, n1, n2 = p
     return db = k2 * hil(1, a, n2) - k4 * b
 end
 
@@ -42,8 +40,8 @@ Asymmetric parameter set
 ps1 = (k1=20., k2=20., k3=5., k4=5., n1=4., n2=1.)
 tend = 4.0
 
-sol1 = solve(ODEProblem(model47!, LVector(a=3., b=1.), tend, ps1))
-sol2 = solve(ODEProblem(model47!, LVector(a=1., b=3.), tend, ps1))
+sol1 = solve(ODEProblem(model47!, [3., 1.], tend, ps1))
+sol2 = solve(ODEProblem(model47!, [1., 3.], tend, ps1))
 
 ax1 = plot(sol1, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.7A (1)")
 ax2 = plot(sol2, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.7A (2)")
@@ -93,8 +91,8 @@ Symmetric parameter set
 ps2 = LVector(k1=20., k2=20., k3=5., k4=5., n1=4., n2=4.)
 
 tend = 4.0
-sol1 = solve(ODEProblem(model47!, LVector(a=3., b=1.), tend, ps2))
-sol2 = solve(ODEProblem(model47!, LVector(a=1., b=3.), tend, ps2))
+sol1 = solve(ODEProblem(model47!, [3., 1.], tend, ps2))
+sol2 = solve(ODEProblem(model47!, [1., 3.], tend, ps2))
 
 ax1 = plot(sol1, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.8A (1)")
 ax2 = plot(sol2, xlabel="Time", ylabel="Concentration", legend=:right, title= "Fig 4.8A (2)")
diff --git a/docs/figs/ch4-15.jl b/docs/figs/ch4-15.jl
index ef01c490..a7fda06f 100644
--- a/docs/figs/ch4-15.jl
+++ b/docs/figs/ch4-15.jl
@@ -7,8 +7,6 @@ Oscillatory networks.
 ===#
 
 using DifferentialEquations
-using LabelledArrays
-using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
@@ -19,13 +17,13 @@ PNG(fig) = display("image/png", fig)
 
 function dA415(u, p, t)
     a, b = u
-    @unpack k0, k1, k2, n = p
+    k0, k1, k2, n = p
     return dA = k0 - k1 * a * (1 + b^n)
 end
 
 function dB415(u, p, t)
     a, b = u
-    @unpack k0, k1, k2, n = p
+    k0, k1, k2, n = p
     return dB = k1 * a * (1 + b^n) - k2 * b
 end
 
@@ -38,10 +36,10 @@ end
 #---
 ps1 = (k0 = 8., k1 = 1., k2 = 5., n = 2.)
 u0s = (
-    LVector(a=1.5, b=1.0),
-    LVector(a=0.0, b=1.0),
-    LVector(a=0.0, b=3.0),
-    LVector(a=2.0, b=0.0),
+    [1.5, 1.0],
+    [0.0, 1.0],
+    [0.0, 3.0],
+    [2.0, 0.0],
 )
 
 tend = 8.
@@ -98,10 +96,10 @@ Oscillatory parameter set
 ps2 = (k0 = 8., k1 = 1., k2 = 5., n = 2.5)
 tend = 100.0
 u0s = (
-    LVector(a=1.5, b=1.0),
-    LVector(a=0.0, b=1.0),
-    LVector(a=0.0, b=3.0),
-    LVector(a=2.0, b=0.0),
+    [1.5, 1.0],
+    [0.0, 1.0],
+    [0.0, 3.0],
+    [2.0, 0.0],
 )
 
 sols = map(u0s) do u0
@@ -148,7 +146,7 @@ fig |> PNG
 
 # ## Fig 4.17
 
-sol = solve(ODEProblem(model415!, LVector(a=2.0, b=1.5), 10.0, ps2))
+sol = solve(ODEProblem(model415!, [2.0, 1.5], 10.0, ps2))
 
 fig = plot(title="Fig 4.17")
 plot!(fig, sol, idxs=(1, 2), label=nothing, arrow=:closed)
diff --git a/docs/figs/ch4-18.jl b/docs/figs/ch4-18.jl
index b9ecdbde..968fc7c3 100644
--- a/docs/figs/ch4-18.jl
+++ b/docs/figs/ch4-18.jl
@@ -7,8 +7,6 @@ See also [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl
 ===#
 
 using DifferentialEquations
-using SimpleUnPack
-using LabelledArrays
 using Plots
 Plots.default(linewidth=2)
 
@@ -21,15 +19,15 @@ hil(x, k, n) = hil(x^n, k^n)
 
 #---
 function model418!(D, u, p, t)
-    @unpack a, b = u
-    @unpack k1, k2, k3, k4, k5, n = p
+    a, b = u
+    k1, k2, k3, k4, k5, n = p
     D.a = k1 * hil(1, b, n) - (k3 + k5) * a
     D.b = k2 + k5 * a - k4 * b
 end
 
 function ainf(k1)
-    ps = LVector(k1 = k1, k2 = 5., k3 = 5., k4 = 5., k5 = 2., n = 4.)
-    u0 = LVector(a=0., b=0.)
+    ps = (k1 = k1, k2 = 5., k3 = 5., k4 = 5., k5 = 2., n = 4.)
+    u0 = [0., 0.]
     prob = SteadyStateProblem(model418!, u0, ps)
     sol = solve(prob, DynamicSS(Rodas5()))
     return sol.u.a
diff --git a/docs/figs/ch5-10.jl b/docs/figs/ch5-10.jl
index e0a8dce6..0fc5c531 100644
--- a/docs/figs/ch5-10.jl
+++ b/docs/figs/ch5-10.jl
@@ -5,7 +5,6 @@ Methionine model
 ===#
 
 using DifferentialEquations
-using LabelledArrays
 using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
@@ -58,7 +57,7 @@ ps = (
 )
 
 tend = 5.
-u0 = LVector(AdoMet=10., AdoHcy=10.)
+u0 = [AdoMet=10., AdoHcy=10.]
 
 prob = ODEProblem(metmodel!, u0, tend, ps)
 sol = solve(prob)
diff --git a/docs/figs/ch7-11.jl b/docs/figs/ch7-11.jl
index 02f352df..6fda718b 100644
--- a/docs/figs/ch7-11.jl
+++ b/docs/figs/ch7-11.jl
@@ -5,8 +5,6 @@ Model of phage lambda decision switch
 ===#
 
 using DifferentialEquations
-using LabelledArrays
-using SimpleUnPack
 using Plots
 Plots.default(linewidth=2)
 
@@ -19,7 +17,7 @@ function model711(u, p, t)
     r, c = u
     rd = r / 2 ## r Dimer
     cd = c / 2 ## c Dimer
-    @unpack K1, K2, K3, K4 , delta_r, delta_c, a, b = p
+    K1, K2, K3, K4 , delta_r, delta_c, a, b = p
     f1 = K1 * rd^2
     f2 = K2 * rd
     f3 = K3 * cd