Hi everybody,

I am trying to solve the following system of ODE with DifferentialEquations.jl

```
function newton!(t, y, dy, μ)
r = norm(y[1:3])
dy[1:3] = y[4:6]
dy[4:6] = -μ * y[1:3] / r^3
end
```

Without units this works perfectly:

```
function test()
r0 = [1131.340, -2282.343, 6672.423]
v0 = [-5.64305, 4.30333, 2.42879]
Δt = 86400.0*365
mu = 398600.4418
rv0 = [r0; v0]
prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu, rv0, (0.0, Δt))
solve(prob, Vern8())[end]
end
```

The unitful version on the other hand does not:

```
function test2()
r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
mu = 398600.4418u"km^3/s^2"
rv0 = [r0; v0]
prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))
solve(prob, Vern8())[end]
end
```

```
TypeError: setindex!: in typeassert, expected Quantity{Float64, Dimensions:{}, Units:{}}, got Float64
```

The state vector `rv0`

with heterogenous units seems to be the problem.

I realise that I could avoid this problem by defining the system of ODE as separate functions in the non-mutating form but this is not an option because the above will be just a small part of the RHS in the complete system.

Is there a way I can have my cake and also eat it?