# How can I simplify this code

**URL:** https://discourse.julialang.org/t/how-can-i-simplify-this-code/107757
**Category:** Modelling & Simulations
**Tags:** question, modelingtoolkit
**Created:** [December 18, 2023, 8:02am UTC](https://discourse.julialang.org/t/how-can-i-simplify-this-code/107757 "2023-12-18T08:02:18Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 18, 2023, 8:02am UTC](https://discourse.julialang.org/t/how-can-i-simplify-this-code/107757/1 "2023-12-18T08:02:18Z")

</div>

How can I simplify the creation of a system of equations in a loop:

```julia
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra

G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration [m/s²]
L0::Float64 = 5.0 # initial segment length [m]
V0::Float64 = 2 # initial velocity of lowest mass [m/s]
M0::Float64 = 0.5 # mass per particle [kg]
C_SPRING::Float64 = 50 # spring constant
segments::Int64 = 5 # number of tether segments [-]
α0 = π/10 # initial tether angle [rad]
duration = 30.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
    v0 = (i-1)*V0/segments
    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 t 
@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] = zeros(3, segments)
D = Differential(t)

eqs1 = vcat(D.(pos) ~ vel,
            D.(vel) ~ acc)
eqs2 = []
for i in segments:-1:1
    global eqs2
    eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
    eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
    eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
    eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
    eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
    eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
    eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
    if i == segments
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
    else
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
    end
    eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
eqs2 = vcat(eqs2, acc[:, 1] .~ zeros(3))
eqs = vcat(eqs1..., eqs2)
     
@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

```

For example:  
Can I use only one vcat call for these two lines?  
Or can I use some other function than vcat?

```julia
    eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
    eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])

```

I tried `push!` and `append!` with no success…

I found out, some of the equations in the array eqs2 have the type `Equation` and some the type `BasicSymbol`. Perhaps that causes the issues?

Why is this not of type Equation:

```julia
julia> eqs2[1]
(broadcast(~, (segment(t))[Colon(), 5], broadcast(-, (pos(t))[Colon(), 6], (pos(t))[Colon(), 5])))[1]

julia> typeof(eqs2[1])
SymbolicUtils.BasicSymbolic{Any}

```

But this is:

```julia
julia> eqs2[end]
damping(t) ~ 473.0 / ((1//5)*length(t))

julia> typeof(eqs2[end])
Equation

```