I have two coupled, first order, non-linear ODEs that I’m trying to solve. However, both of the two equations contain additional derivatives on the right-hand sides of the equation (and there don’t seem to be any obvious algebraic means to rectify that). DifferentialEquations.jl requires that the equations be in the form,

\frac{du}{dt} = f(u, p, t) ,

so I’m not sure I can use DifferenetialEquations.jl to solve my ODE as it is presently formulated. I found one example problem that has derivatives on the right side: the three-body problem. Copying from the DiffEq docs, the three body problem is formulated as follows:

\frac{dy_1}{dt} = y_1 + 2\frac{dy_2}{dt} - \bar{μ}\frac{y_1+μ}{D_1} - μ\frac{y_1-\bar{μ}}{D_2}

\frac{dy_2}{dt} = y_2 - 2\frac{dy_1}{dt} - \bar{μ}\frac{y_2}{D_1} - μ\frac{y_2}{D_2}

However, the implementation that I see on Github for the right-hand side function `f(u, p, t)`

doesn’t make much sense to me. Here is the implementation on Github:

```
const threebody_μ = big(0.012277471); const threebody_μ′ = 1 - threebody_μ
threebody = (du,u,p,t) -> begin
# 1 = y₁
# 2 = y₂
# 3 = y₁'
# 4 = y₂'
D₁ = ((u[1]+threebody_μ)^2 + u[2]^2)^(3/2)
D₂ = ((u[1]-threebody_μ′)^2 + u[2]^2)^(3/2)
du[1] = u[3]
du[2] = u[4]
du[3] = u[1] + 2u[4] - threebody_μ′*(u[1]+threebody_μ)/D₁ - threebody_μ*(u[1]-threebody_μ′)/D₂
du[4] = u[2] - 2u[3] - threebody_μ′*u[2]/D₁ - threebody_μ*u[2]/D₂
end
```

It appears to me that the lines defining `du[3]`

and `du[4]`

are actually defining the second derivatives of y_1 and y_2, rather than the first derivatives that appear in the ODE.

So, I have a few questions:

- Is this implementation correct?
- If it is correct, can anyone explain to me why it is correct and how it works?
- If it’s incorrect, is there a correct way that one can do this sort of thing with DifferentialEquations.jl?

Thanks in advance.