diff --git a/Project.toml b/Project.toml
index 4f1ea3d..fe7ad03 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ Conda = "1.10.2"
 CondaPkg = "0.2.23"
 ControlPlots = "0.2.1"
 LiveServer = "1.3.1"
-ModelingToolkit = "~9.38.0"
+ModelingToolkit = "~9.40.1"
 OrdinaryDiffEq = "6.87.0"
 PackageCompiler = "2.1.17"
 Parameters = "0.12.3"
diff --git a/src/Tether_01.jl b/src/Tether_01.jl
index 972623a..11ea634 100644
--- a/src/Tether_01.jl
+++ b/src/Tether_01.jl
@@ -1,5 +1,5 @@
 # Example one: Falling mass.
-using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit, OrdinaryDiffEq, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -7,7 +7,7 @@ G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration
 # definiting the model
 @variables   pos(t)[1:3] = [0.0, 0.0,  0.0]
 @variables   vel(t)[1:3] = [0.0, 0.0, 50.0] 
-@variables   acc(t)[1:3] = [0.0, 0.0, -9.81] 
+@variables   acc(t)[1:3]
 
 eqs = vcat(D(pos) ~ vel,
            D(vel) ~ acc,
diff --git a/src/Tether_02.jl b/src/Tether_02.jl
index 0b66f7b..b715ff8 100644
--- a/src/Tether_02.jl
+++ b/src/Tether_02.jl
@@ -1,5 +1,5 @@
 # Example two: Falling mass, attached to linear spring
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -9,10 +9,10 @@ L0::Float64 = -10.0                             # initial spring length      [m]
 @parameters mass=1.0 c_spring=50.0 damping=0.5 l0=L0
 @variables   pos(t)[1:3] = [0.0, 0.0,  L0]
 @variables   vel(t)[1:3] = [0.0, 0.0,  0.0] 
-@variables   acc(t)[1:3] = G_EARTH
-@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
-@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
-@variables norm1(t) = abs(l0) spring_vel(t) = 0.0
+@variables   acc(t)[1:3]
+@variables unit_vector(t)[1:3]
+@variables spring_force(t)[1:3]
+@variables norm1(t) spring_vel(t)
 
 eqs = vcat(D(pos)       ~ vel,
            D(vel)       ~ acc,
diff --git a/src/Tether_03.jl b/src/Tether_03.jl
index a6da4a6..dc4c566 100644
--- a/src/Tether_03.jl
+++ b/src/Tether_03.jl
@@ -1,6 +1,6 @@
 # Example three: Falling mass, attached to non-linear spring without compression stiffness,
 # initially moving upwards with 4 m/s.
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -11,11 +11,11 @@ V0::Float64 = 4                                 # initial velocity           [m/
 @parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
 @variables   pos(t)[1:3] = [0.0, 0.0,  L0]
 @variables   vel(t)[1:3] = [0.0, 0.0,  V0] 
-@variables   acc(t)[1:3] = G_EARTH
-@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
-@variables c_spring(t) = c_spring0
-@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
-@variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
+@variables   acc(t)[1:3]
+@variables unit_vector(t)[1:3]
+@variables c_spring(t)
+@variables spring_force(t)[1:3]
+@variables force(t) norm1(t) spring_vel(t)
 
 eqs = vcat(D(pos)       ~ vel,
            D(vel)       ~ acc,
diff --git a/src/Tether_03b.jl b/src/Tether_03b.jl
index f862b0b..ae3b128 100644
--- a/src/Tether_03b.jl
+++ b/src/Tether_03b.jl
@@ -1,6 +1,6 @@
 # Example three: Falling mass, attached to non-linear spring without compression stiffness
 # with callback, initially moving upwards with 4 m/s.
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -11,11 +11,11 @@ V0::Float64 = 4                                 # initial velocity           [m/
 @parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
 @variables   pos(t)[1:3] = [0.0, 0.0,  L0]
 @variables   vel(t)[1:3] = [0.0, 0.0,  V0] 
-@variables   acc(t)[1:3] = G_EARTH
-@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
-@variables c_spring(t) = c_spring0
-@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
-@variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
+@variables   acc(t)[1:3]
+@variables unit_vector(t)[1:3]
+@variables c_spring(t)
+@variables spring_force(t)[1:3]
+@variables force(t) norm1(t) spring_vel(t)
 
 eqs = vcat(D.(pos)      ~ vel,
            D.(vel)      ~ acc,
diff --git a/src/Tether_03c.jl b/src/Tether_03c.jl
index ff0edae..d4c95c8 100644
--- a/src/Tether_03c.jl
+++ b/src/Tether_03c.jl
@@ -1,6 +1,6 @@
 # Example three: Falling mass, attached to non-linear spring without compression stiffness,
 # initially moving upwards with 4 m/s. Comparing results with and without callbacks.
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -12,11 +12,11 @@ function model3(G_EARTH, L0, V0)
     @parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
     @variables   pos(t)[1:3] = [0.0, 0.0,  L0]
     @variables   vel(t)[1:3] = [0.0, 0.0,  V0] 
-    @variables   acc(t)[1:3] = G_EARTH
-    @variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
-    @variables c_spring(t) = c_spring0
-    @variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
-    @variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
+    @variables   acc(t)[1:3]
+    @variables unit_vector(t)[1:3]
+    @variables c_spring(t)
+    @variables spring_force(t)[1:3]
+    @variables force(t) norm1(t) spring_vel(t)
 
     eqs = vcat(D.(pos)   ~ vel,
             D.(vel)      ~ acc,
diff --git a/src/Tether_04.jl b/src/Tether_04.jl
index 79a10ff..e91e462 100644
--- a/src/Tether_04.jl
+++ b/src/Tether_04.jl
@@ -5,7 +5,7 @@ for l < l_0) and n tether segments.
 # TODO: Distribute force correctly
 # TODO: Add 2D plot
 
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -14,33 +14,23 @@ V0::Float64 = 4                                 # initial velocity of lowest mas
 segments::Int64 = 2                             # number of tether segments         [-]
 POS0 = zeros(3, segments+1)
 VEL0 = zeros(3, segments+1)
-ACC0 = zeros(3, segments+1)
-SEGMENTS0 = zeros(3, segments) 
-UNIT_VECTORS0 = zeros(3, segments)
 for i in 1:segments+1
     POS0[:, i] .= [0.0, 0, -(i-1)*L0]
     VEL0[:, i] .= [0.0, 0, (i-1)*V0/segments]
 end
-for i in 2:segments+1
-    ACC0[:, i] .= G_EARTH
-end
-for i in 1:segments
-    UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-    SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-end
 
 # defining the model, Z component upwards
 @parameters mass=1.0 c_spring0=50.0 damping=0.5 l_seg=L0
 @variables pos(t)[1:3, 1:segments+1]  = POS0
 @variables vel(t)[1:3, 1:segments+1]  = VEL0
-@variables acc(t)[1:3, 1:segments+1]  = ACC0
-@variables segment(t)[1:3, 1:segments]  = SEGMENTS0
-@variables unit_vector(t)[1:3, 1:segments]  = UNIT_VECTORS0
-@variables norm1(t)[1:segments] = l_seg * ones(segments)
-@variables rel_vel(t)[1:3, 1:segments]  = zeros(3, segments)
-@variables spring_vel(t)[1:segments] = zeros(segments)
-@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
-@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
+@variables acc(t)[1:3, 1:segments+1]
+@variables segment(t)[1:3, 1:segments]
+@variables unit_vector(t)[1:3, 1:segments]
+@variables norm1(t)[1:segments]
+@variables rel_vel(t)[1:3, 1:segments]
+@variables spring_vel(t)[1:segments]
+@variables c_spring(t)[1:segments]
+@variables spring_force(t)[1:3, 1:segments]
 
 # basic differential equations
 eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_05.jl b/src/Tether_05.jl
index 926a79f..bd720e5 100644
--- a/src/Tether_05.jl
+++ b/src/Tether_05.jl
@@ -1,6 +1,6 @@
 # Tutorial example simulating a 3D mass-spring system with a nonlinear spring (no spring forces
 # for l < l_0) and n tether segments. 
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/s²]
@@ -13,9 +13,6 @@ segments::Int64 = 5                             # number of tether segments
 duration = 10.0                                 # duration of the simulation        [s]
 POS0 = zeros(3, segments+1)
 VEL0 = zeros(3, segments+1)
-ACC0 = zeros(3, segments+1)
-SEGMENTS0 = zeros(3, segments) 
-UNIT_VECTORS0 = zeros(3, segments)
 for i in 1:segments+1
     local l0
     l0 = -(i-1) * L0
@@ -23,27 +20,20 @@ for i in 1:segments+1
     POS0[:, i] .= [sin(α0) * l0, 0, cos(α0) * l0]
     VEL0[:, i] .= [sin(α0) * v0, 0, cos(α0) * v0]
 end
-for i in 2:segments+1
-    ACC0[:, i] .= G_EARTH
-end
-for i in 1:segments
-    UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-    SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-end
 
 # defining the model, Z component upwards
 @parameters mass=M0 c_spring0=C_SPRING damping=0.5 l_seg=L0
 @variables pos(t)[1:3, 1:segments+1]  = POS0
 @variables vel(t)[1:3, 1:segments+1]  = VEL0
-@variables acc(t)[1:3, 1:segments+1]  = ACC0
-@variables segment(t)[1:3, 1:segments]  = SEGMENTS0
-@variables unit_vector(t)[1:3, 1:segments]  = UNIT_VECTORS0
-@variables norm1(t)[1:segments] = l_seg * ones(segments)
-@variables rel_vel(t)[1:3, 1:segments]  = zeros(3, segments)
-@variables spring_vel(t)[1:segments] = zeros(segments)
-@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
-@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
-@variables total_force(t)[1:3, 1:segments+1] = zeros(3, segments+1)
+@variables acc(t)[1:3, 1:segments+1]
+@variables segment(t)[1:3, 1:segments]
+@variables unit_vector(t)[1:3, 1:segments]
+@variables norm1(t)[1:segments]
+@variables rel_vel(t)[1:3, 1:segments]
+@variables spring_vel(t)[1:segments]
+@variables c_spring(t)[1:segments]
+@variables spring_force(t)[1:3, 1:segments]
+@variables total_force(t)[1:3, 1:segments+1]
 
 # basic differential equations
 eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_06.jl b/src/Tether_06.jl
index a2207af..efacbc6 100644
--- a/src/Tether_06.jl
+++ b/src/Tether_06.jl
@@ -1,7 +1,7 @@
 # Tutorial example simulating a 3D mass-spring system with a nonlinear spring (1% spring forces
 # for l < l_0), n tether segments and reel-in and reel-out. The compression stiffness is hardcoded
 # to 1% of the spring stiffness.
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ControlPlots
 
@@ -21,37 +21,29 @@ mass_per_meter::Float64 = RHO_TETHER * π * (D_TETHER/2000.0)^2
 # calculating consistant initial conditions
 POS0 = zeros(3, SEGMENTS+1)
 VEL0 = zeros(3, SEGMENTS+1)
-ACC0 = zeros(3, SEGMENTS+1)
-SEGMENTS0 = zeros(3, SEGMENTS) 
-UNIT_VECTORS0 = zeros(3, SEGMENTS)
 for i in 1:SEGMENTS+1
     l0_ = -(i-1)*L0/SEGMENTS
     POS0[:, i] .= [sin(α0) * l0_, 0, cos(α0) * l0_]
     VEL0[:, i] .= [0.0, 0, 0]
 end
-for i in 1:SEGMENTS
-    ACC0[:, i+1] .= G_EARTH
-    UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-    SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-end
 
 # defining the model, Z component upwards
 @parameters c_spring0=C_SPRING/(L0/SEGMENTS) l_seg=L0/SEGMENTS
 @variables pos(t)[1:3, 1:SEGMENTS+1]  = POS0
 @variables vel(t)[1:3, 1:SEGMENTS+1]  = VEL0
-@variables acc(t)[1:3, 1:SEGMENTS+1]  = ACC0
-@variables segment(t)[1:3, 1:SEGMENTS]  = SEGMENTS0
-@variables unit_vector(t)[1:3, 1:SEGMENTS]  = UNIT_VECTORS0
-@variables len(t) = L0
-@variables c_spring(t) = c_spring0
-@variables damping(t) = DAMPING  / l_seg
-@variables m_tether_particle(t) = mass_per_meter * l_seg
-@variables norm1(t)[1:SEGMENTS] = l_seg * ones(SEGMENTS)
-@variables rel_vel(t)[1:3, 1:SEGMENTS]  = zeros(3, SEGMENTS)
-@variables spring_vel(t)[1:SEGMENTS] = zeros(SEGMENTS)
-@variables c_spr(t)[1:SEGMENTS] = c_spring0 * ones(SEGMENTS)
-@variables spring_force(t)[1:3, 1:SEGMENTS] = zeros(3, SEGMENTS)
-@variables total_force(t)[1:3, 1:SEGMENTS+1] = zeros(3, SEGMENTS+1)
+@variables acc(t)[1:3, 1:SEGMENTS+1]
+@variables segment(t)[1:3, 1:SEGMENTS]
+@variables unit_vector(t)[1:3, 1:SEGMENTS]
+@variables len(t)
+@variables c_spring(t)
+@variables damping(t)
+@variables m_tether_particle(t)
+@variables norm1(t)[1:SEGMENTS]
+@variables rel_vel(t)[1:3, 1:SEGMENTS]
+@variables spring_vel(t)[1:SEGMENTS]
+@variables c_spr(t)[1:SEGMENTS]
+@variables spring_force(t)[1:3, 1:SEGMENTS]
+@variables total_force(t)[1:3, 1:SEGMENTS+1]
 
 # basic differential equations
 eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_06b.jl b/src/Tether_06b.jl
index f370962..4e4c281 100644
--- a/src/Tether_06b.jl
+++ b/src/Tether_06b.jl
@@ -1,7 +1,7 @@
 # Tutorial example simulating a 3D mass-spring system with a nonlinear spring (no spring forces
 # for l < l_0), n tether segments and reel-in and reel-out. Can create a video of the simulation.
 # For creating the video, set save=true in the Settings struct.
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ControlPlots
 
@@ -22,41 +22,33 @@ end
 function calc_initial_state(se)
     POS0 = zeros(3, se.segments+1)
     VEL0 = zeros(3, se.segments+1)
-    ACC0 = zeros(3, se.segments+1)
-    SEGMENTS0 = zeros(3, se.segments) 
-    UNIT_VECTORS0 = zeros(3, se.segments)
     for i in 1:se.segments+1
         l0 = -(i-1)*se.l0/se.segments
         POS0[:, i] .= [sin(se.α0) * l0, 0, cos(se.α0) * l0]
         VEL0[:, i] .= [0, 0, 0]
     end
-    for i in 1:se.segments
-        ACC0[:, i+1] .= se.g_earth
-        UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-        SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-    end
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0
+    POS0, VEL0
 end
 
 function model(se)
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0 = calc_initial_state(se)
+    POS0, VEL0 = calc_initial_state(se)
     mass_per_meter = se.rho_tether * π * (se.d_tether/2000.0)^2
     @parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
     @variables pos(t)[1:3, 1:se.segments+1]  = POS0
     @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
-    @variables acc(t)[1:3, 1:se.segments+1]  = ACC0
-    @variables segment(t)[1:3, 1:se.segments]  = SEGMENTS0
-    @variables unit_vector(t)[1:3, 1:se.segments]  = UNIT_VECTORS0
-    @variables len(t) = se.l0
-    @variables c_spring(t) = c_spring0
-    @variables damping(t) = se.damping  / l_seg
-    @variables m_tether_particle(t) = mass_per_meter * l_seg
-    @variables norm1(t)[1:se.segments] = l_seg * ones(se.segments)
-    @variables rel_vel(t)[1:3, 1:se.segments]  = zeros(3, se.segments)
-    @variables spring_vel(t)[1:se.segments] = zeros(se.segments)
-    @variables c_spr(t)[1:se.segments] = c_spring0 * ones(se.segments)
-    @variables spring_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables total_force(t)[1:3, 1:se.segments+1] = zeros(3, se.segments+1)
+    @variables acc(t)[1:3, 1:se.segments+1]
+    @variables segment(t)[1:3, 1:se.segments]
+    @variables unit_vector(t)[1:3, 1:se.segments]
+    @variables len(t)
+    @variables c_spring(t)
+    @variables damping(t)
+    @variables m_tether_particle(t)
+    @variables norm1(t)[1:se.segments]
+    @variables rel_vel(t)[1:3, 1:se.segments]
+    @variables spring_vel(t)[1:se.segments]
+    @variables c_spr(t)[1:se.segments]
+    @variables spring_force(t)[1:3, 1:se.segments]
+    @variables total_force(t)[1:3, 1:se.segments+1]
 
     # basic differential equations
     eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_06c.jl b/src/Tether_06c.jl
index 010b4b5..1360a5b 100644
--- a/src/Tether_06c.jl
+++ b/src/Tether_06c.jl
@@ -1,6 +1,6 @@
 # Tutorial example simulating a 3D mass-spring system with a nonlinear spring (no spring forces
 # for l < l_0), n tether segments, reel-in and reel-out and continues callbacks. 
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ControlPlots
 
@@ -25,43 +25,35 @@ end
 function calc_initial_state(se)
     POS0 = zeros(3, se.segments+1)
     VEL0 = zeros(3, se.segments+1)
-    ACC0 = zeros(3, se.segments+1)
-    SEGMENTS0 = zeros(3, se.segments) 
-    UNIT_VECTORS0 = zeros(3, se.segments)
     for i in 1:se.segments+1
         l0 = -(i-1)*se.l0/se.segments
         v0 = (i-1)*se.v0/se.segments
         POS0[:, i] .= [sin(se.α0) * l0, 0, cos(se.α0) * l0]
         VEL0[:, i] .= [sin(se.α0) * v0, 0, cos(se.α0) * v0]
     end
-    for i in 1:se.segments
-        ACC0[:, i+1] .= se.g_earth
-        UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-        SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-    end
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0
+    POS0, VEL0
 end
 
 function model(se)
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0 = calc_initial_state(se)
+    POS0, VEL0 = calc_initial_state(se)
     mass_per_meter = se.rho_tether * π * (se.d_tether/2000.0)^2
     @parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
     @parameters rel_compression_stiffness = se.rel_compression_stiffness    
     @variables pos(t)[1:3, 1:se.segments+1]  = POS0
     @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
-    @variables acc(t)[1:3, 1:se.segments+1]  = ACC0
-    @variables segment(t)[1:3, 1:se.segments]  = SEGMENTS0
-    @variables unit_vector(t)[1:3, 1:se.segments]  = UNIT_VECTORS0
-    @variables length(t) = se.l0
-    @variables c_spring(t) = c_spring0
-    @variables damping(t) = se.damping  / l_seg
-    @variables m_tether_particle(t) = mass_per_meter * l_seg
-    @variables norm1(t)[1:se.segments] = l_seg * ones(se.segments)
-    @variables rel_vel(t)[1:3, 1:se.segments]  = zeros(3, se.segments)
-    @variables spring_vel(t)[1:se.segments] = zeros(se.segments)
-    @variables c_spr(t)[1:se.segments] = c_spring0 * ones(se.segments)
-    @variables spring_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables total_force(t)[1:3, 1:se.segments+1] = zeros(3, se.segments+1)
+    @variables acc(t)[1:3, 1:se.segments+1]
+    @variables segment(t)[1:3, 1:se.segments]
+    @variables unit_vector(t)[1:3, 1:se.segments]
+    @variables length(t)
+    @variables c_spring(t)
+    @variables damping(t)
+    @variables m_tether_particle(t)
+    @variables norm1(t)[1:se.segments]
+    @variables rel_vel(t)[1:3, 1:se.segments]
+    @variables spring_vel(t)[1:se.segments]
+    @variables c_spr(t)[1:se.segments]
+    @variables spring_force(t)[1:3, 1:se.segments]
+    @variables total_force(t)[1:3, 1:se.segments+1]
 
     # basic differential equations
     eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_07.jl b/src/Tether_07.jl
index 90a29aa..14115b2 100644
--- a/src/Tether_07.jl
+++ b/src/Tether_07.jl
@@ -1,6 +1,6 @@
 # Tutorial example simulating a 3D mass-spring system with a nonlinear spring (1% stiffness
 # for l < l_0), n tether segments, tether drag and reel-in and reel-out. 
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ControlPlots
 
@@ -31,47 +31,39 @@ end
 function calc_initial_state(se)
     POS0 = zeros(3, se.segments+1)
     VEL0 = zeros(3, se.segments+1)
-    ACC0 = zeros(3, se.segments+1)
-    SEGMENTS0 = zeros(3, se.segments) 
-    UNIT_VECTORS0 = zeros(3, se.segments)
     for i in 1:se.segments+1
         l0 = -(i-1)*se.l0/se.segments
         v0 = -(i-1)*se.v_ro/se.segments
         POS0[:, i] .= [sin(se.α0) * l0, 0, cos(se.α0) * l0]
         VEL0[:, i] .= [sin(se.α0) * v0, 0, cos(se.α0) * v0]
     end
-    for i in 1:se.segments
-        ACC0[:, i+1] .= se.g_earth
-        UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
-        SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-    end
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0
+    POS0, VEL0
 end
 
 function model(se)
-    POS0, VEL0, ACC0, SEGMENTS0, UNIT_VECTORS0 = calc_initial_state(se)
+    POS0, VEL0 = calc_initial_state(se)
     mass_per_meter = se.rho_tether * π * (se.d_tether/2000.0)^2
     @parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
     @parameters rel_compression_stiffness = se.rel_compression_stiffness
     @variables pos(t)[1:3, 1:se.segments+1]  = POS0
     @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
-    @variables acc(t)[1:3, 1:se.segments+1]  = ACC0
-    @variables segment(t)[1:3, 1:se.segments]  = SEGMENTS0
-    @variables unit_vector(t)[1:3, 1:se.segments]  = UNIT_VECTORS0
-    @variables length(t) = se.l0
-    @variables c_spring(t) = c_spring0
-    @variables damping(t) = se.damping  / l_seg
-    @variables m_tether_particle(t) = mass_per_meter * l_seg
-    @variables norm1(t)[1:se.segments] = l_seg * ones(se.segments)
-    @variables rel_vel(t)[1:3, 1:se.segments]  = zeros(3, se.segments)
-    @variables spring_vel(t)[1:se.segments] = zeros(se.segments)
-    @variables c_spr(t)[1:se.segments] = c_spring0 * ones(se.segments)
-    @variables spring_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables v_apparent(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables v_app_perp(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables norm_v_app(t)[1:se.segments] = ones(se.segments)
-    @variables half_drag_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables total_force(t)[1:3, 1:se.segments+1] = zeros(3, se.segments+1)
+    @variables acc(t)[1:3, 1:se.segments+1]
+    @variables segment(t)[1:3, 1:se.segments]
+    @variables unit_vector(t)[1:3, 1:se.segments]
+    @variables length(t)
+    @variables c_spring(t)
+    @variables damping(t)
+    @variables m_tether_particle(t)
+    @variables norm1(t)[1:se.segments]
+    @variables rel_vel(t)[1:3, 1:se.segments]
+    @variables spring_vel(t)[1:se.segments]
+    @variables c_spr(t)[1:se.segments]
+    @variables spring_force(t)[1:3, 1:se.segments]
+    @variables v_apparent(t)[1:3, 1:se.segments]
+    @variables v_app_perp(t)[1:3, 1:se.segments]
+    @variables norm_v_app(t)[1:se.segments]
+    @variables half_drag_force(t)[1:3, 1:se.segments]
+    @variables total_force(t)[1:3, 1:se.segments+1]
 
     # basic differential equations
     eqs1 = vcat(D.(pos) .~ vel,
diff --git a/src/Tether_08.jl b/src/Tether_08.jl
index 8e19fc3..bd97221 100644
--- a/src/Tether_08.jl
+++ b/src/Tether_08.jl
@@ -2,7 +2,7 @@
 # for l < l_0), n tether segments, tether drag and reel-in and reel-out. 
 # New feature: A steady state solver is used to find the initial tether shape for any
 # given pair of endpoints, which is then used as the initial condition for the simulation.
-using ModelingToolkit, OrdinaryDiffEq, SteadyStateDiffEq, LinearAlgebra, Timers, Parameters
+using ModelingToolkit, OrdinaryDiffEq, SteadyStateDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots
 tic()
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using ControlPlots
@@ -41,16 +41,12 @@ function calc_initial_state(se; p1, p2)
     end
     POS0 = zeros(3, se.segments+1)
     VEL0 = zeros(3, se.segments+1)
-    ACC0 = zeros(3, se.segments+1)
     # use a linear interpolation between p1 and p2 for the intermediate points
     for i in 1:se.segments+1
         Δ = (p2-p1) / se.segments
         POS0[:, i] .= p1 + (i-1) * Δ
     end
-    for i in 1:se.segments
-        ACC0[:, i+1] .= se.g_earth
-    end
-    POS0, VEL0, ACC0
+    POS0, VEL0
 end
 
 function model(se; p1=[0,0,0], p2=nothing, fix_p1=true, fix_p2=false)
@@ -68,11 +64,11 @@ function model(se; p1=[0,0,0], p2=nothing, fix_p1=true, fix_p2=false)
         @assert ! fix_p2                || error("if p2 undefined it cannot be fixed")
     end
     # straight line approximation for the tether
-    POS0, VEL0, ACC0 = calc_initial_state(se; p1, p2)
+    POS0, VEL0 = calc_initial_state(se; p1, p2)
     # find steady state
     v_ro = se.v_ro      # save the reel-out speed
     se.v_ro = 0         # v_ro must be zero, otherwise finding the steady state is not possible
-    simple_sys, pos, =  model(se, p1, p2, true, true, POS0, VEL0, ACC0)
+    simple_sys, pos, =  model(se, p1, p2, true, true, POS0, VEL0)
     tspan = (0.0, se.duration)
     prob = ODEProblem(simple_sys, nothing, tspan)
     prob1 = SteadyStateProblem(prob)
@@ -80,35 +76,31 @@ function model(se; p1=[0,0,0], p2=nothing, fix_p1=true, fix_p2=false)
     POS0 = sol1[pos]
     # create the real model
     se.v_ro = v_ro
-    model(se, p1, p2, fix_p1, fix_p2, POS0, VEL0, ACC0)
+    model(se, p1, p2, fix_p1, fix_p2, POS0, VEL0)
 end
-function model(se, p1, p2, fix_p1, fix_p2, POS0, VEL0, ACC0)
-    SEGMENTS0 = zeros(3, se.segments)
-    for i in 1:se.segments
-        SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
-    end
+function model(se, p1, p2, fix_p1, fix_p2, POS0, VEL0)
     mass_per_meter = se.rho_tether * π * (se.d_tether/2000.0)^2
     @parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
     @parameters rel_compression_stiffness = se.rel_compression_stiffness
     @variables pos(t)[1:3, 1:se.segments+1]  = POS0
     @variables vel(t)[1:3, 1:se.segments+1]  = VEL0
-    @variables acc(t)[1:3, 1:se.segments+1]  = ACC0
-    @variables segment(t)[1:3, 1:se.segments]  = SEGMENTS0
-    @variables unit_vector(t)[1:3, 1:se.segments]  = zeros(3, se.segments).+[0,0,1]
-    @variables len(t) = se.l0
-    @variables c_spring(t) = c_spring0
-    @variables damping(t) = se.damping  / l_seg
-    @variables m_tether_particle(t) = mass_per_meter * l_seg
-    @variables norm1(t)[1:se.segments] = l_seg * ones(se.segments)
-    @variables rel_vel(t)[1:3, 1:se.segments]  = zeros(3, se.segments)
-    @variables spring_vel(t)[1:se.segments] = zeros(se.segments)
-    @variables c_spr(t)[1:se.segments] = c_spring0 * ones(se.segments)
-    @variables spring_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables v_apparent(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables v_app_perp(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables norm_v_app(t)[1:se.segments] = ones(se.segments)
-    @variables half_drag_force(t)[1:3, 1:se.segments] = zeros(3, se.segments)
-    @variables total_force(t)[1:3, 1:se.segments+1] = zeros(3, se.segments+1)
+    @variables acc(t)[1:3, 1:se.segments+1]
+    @variables segment(t)[1:3, 1:se.segments]
+    @variables unit_vector(t)[1:3, 1:se.segments]
+    @variables len(t)
+    @variables c_spring(t)
+    @variables damping(t)
+    @variables m_tether_particle(t)
+    @variables norm1(t)[1:se.segments]
+    @variables rel_vel(t)[1:3, 1:se.segments]
+    @variables spring_vel(t)[1:se.segments]
+    @variables c_spr(t)[1:se.segments]
+    @variables spring_force(t)[1:3, 1:se.segments]
+    @variables v_apparent(t)[1:3, 1:se.segments]
+    @variables v_app_perp(t)[1:3, 1:se.segments]
+    @variables norm_v_app(t)[1:se.segments]
+    @variables half_drag_force(t)[1:3, 1:se.segments]
+    @variables total_force(t)[1:3, 1:se.segments+1]
 
     # basic differential equations
     eqs1 = vcat(D.(pos) .~ vel,