thank you for the insight to understand this.

after doing some more exploration I found the following code to over come the additional equation.

```
using Symbolics, Plots, DifferentialEquations, ModelingToolkit
n = 2 # number of general coordinates
vars = @variables t (q(t))[1:n] # q[1] = x, q[2] = θ
params = @parameters m₁ m₂ l g F
Q = [F*cos(t), 0]# general forces
t = vars[1]
qs = vars[2]
dt = Differential(t)
dqs = [Differential(q) for q in qs]
dq̇s = [Differential(dt(q)) for q in qs]
q̇s = [dt(q) for q in qs]
v = q̇s[2]*l
T = 0.5*(m₁+m₂)*q̇s[1]^2 + 0.5*m₂*v^2 + m₂*l*q̇s[1]*q̇s[2]*cos(qs[2])
U = -m₂*g*l*cos(qs[2])
L = T - U
Eᵢ = Vector{Num}()
Eᵢᵢ = Vector{Num}()
for (ind, (q, dq̇, dq)) in enumerate(zip(qs, dq̇s, dqs))
push!(Eᵢ, expand_derivatives(dt(expand_derivatives(dq̇(L)))) - expand_derivatives(dq(L)))
push!(Eᵢᵢ, Symbolics.solve_for(Eᵢ[ind] ~ Q[ind], dt(dt(q))))
end
system = [dt(dt(q)) ~ E for (E, q) in zip(Eᵢᵢ,qs)]
tmp = Symbolics.solve_for(system, [dt(dt(qs[1])), dt(dt(qs[2]))])
system = [dt(dt(q)) ~ e for (e, q) in zip(tmp,qs)]
@named HO = ODESystem(system, t, qs, params)
HO = structural_simplify(HO)
X_0 = [qs[1] => 0,
qs[2] => π/6,
q̇s[1] => 0,
q̇s[2] => 0]
M = [ m₁ => 30, # kg
m₂ => 10, # kg
l => 30, # m
g => 9.81, # m/s^2
F => 0] # N
tspan = (0.0, 100.0)
prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)
```

As I understand it `system = [dt(dt(q)) ~ E for (E, q) in zip(Eᵢᵢ,qs)]`

produce the following equations:

```
2-element Vector{Equation}:
Differential(t)(Differential(t)((q(t))[1])) ~ (l*m₂*Differential(t)(Differential(t)((q(t))[2]))*cos((q(t))[2]) - F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂)
Differential(t)(Differential(t)((q(t))[2])) ~ (l*m₂*Differential(t)(Differential(t)((q(t))[1]))*cos((q(t))[2]) + g*l*m₂*sin((q(t))[2])) / (-m₂*(l^2))
```

which find the second derivative of each state as function of the second derivative of the other state

to over come this I add additional step :

```
tmp = Symbolics.solve_for(system, [dt(dt(qs[1])), dt(dt(qs[2]))])
system = [dt(dt(q)) ~ e for (e, q) in zip(tmp,qs)]
```

This step created 2 equations in which the second derivative of each state is a function of lower order derivates of all the other states, which is the “normal” way to write those set of equations.

```
2-element Vector{Equation}:
Differential(t)(Differential(t)((q(t))[1])) ~ (-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂) + (l*m₂*((g*sin((q(t))[2])) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂))*cos((q(t))[2])) / l)*cos((q(t))[2])) / (((-m₂*(cos((q(t))[2])^2)) / (-m₁ - m₂) - 1)*(-m₁ - m₂))
Differential(t)(Differential(t)((q(t))[2])) ~ ((g*sin((q(t))[2])) / l + (((-F*cos(t) - l*m₂*(Differential(t)((q(t))[2])^2)*sin((q(t))[2])) / (-m₁ - m₂))*cos((q(t))[2])) / l) / ((-m₂*(cos((q(t))[2])^2)) / (-m₁ - m₂) - 1)
```

(A lot of text but the important thing is that there aren’t any second derivative at the right hand side of each equation).

Now when I call

```
HO = structural_simplify(HO)
```

I get 4 equations and the following solution.

thank you all for the help.