So, as I mentioned in this topic, I was trying to implement the Lagrangian approach from scratch in Julia, and inquired for existing ways to do that. As the comments on that discussion point out, nothing of the sort exists yet, but the different packages suggested were `Symbolics.jl`

, `ModelingToolkit.jl`

, `DifferentialEquations.jl`

.

I thought about it a bit and came up with a way to use these packages to implement simpler Lagrangians without much hassle.

The general idea is to use the `Symbolics.jl`

to define 3 kinds of derivatives, d/dt, d/dx, d/d\dot{x}. The Lagrangian can be defined in terms of T and V properly, by using `expand_derivatives()`

in conjunction with the above. Then, we can define the Euler - Lagrange equations, use `Symbolics.solve_for()`

to get the double derivatives with which we can define our ODESystem.

An example, for the isotropic Harmonic Oscillator in 2D, is given below:

```
using Symbolics, Plots, DifferentialEquations, ModelingToolkit
@variables t x(t) y(t)
@parameters m, k
d = Differential(t)
D1 = Differential(x)
D2 = Differential(y)
D11 = Differential(d(x))
D22 = Differential(d(y))
T = 0.5*m*(expand_derivatives(d(x))^2 + expand_derivatives(d(y))^2) #KE
V = 0.5*k*(x^2 + y^2) #Potential energy
L = T - V
E₁ = expand_derivatives(d(D11(L))) - expand_derivatives(D1(L))
E₁₁ = Symbolics.solve_for(E₁ ~ 0, d(d(x)))
E₂ = expand_derivatives(d(D22(L))) - expand_derivatives(D2(L))
E₂₂ = Symbolics.solve_for(E₂ ~ 0, d(d(y)))
system = [d(d(x)) ~ E₁₁,
d(d(y)) ~ E₂₂]
@named HO = ODESystem(system, t, [x, y], [m, k])
HO = structural_simplify(HO)
X_0 = [d(x) => 0,
d(y) => 0,
x => 10,
y => 10]
M = [ m => 1
k => 1 ]
tspan = (0.0, 100.0)
prob = ODEProblem(HO, X_0, tspan, M, abstol = 1e-8, reltol = 1e-8)
sol = solve(prob)
```

And, done.

I have plotted this solution below against an expected one.

However, this actually doesn’t solve my original problem, because I have coupled derivatives in there (so the equation for d^{2}x/dt^{2} will contain double derivatives of y and z as well, and same for the other 2 coordinates as well). I would appreciate it if someone can think of a way around that issue.

And, if not, this post will hopefully save someone some time in the future. : )