Following the tutorial, I want to use OrdinaryDiffEq to solve an equation:

using OrdinaryDiffEq
function A(t, p)
ω, ε = p
return [[0.0 1.0]; [-(ω + ε * cos(t)) 0.0]]
end
function mathieu!(du, u, p, t)
du = A(t, p) * u
end
u0 = [1.0, 1.0]
p = [1.0, 1.0]
tspan = (0.0, 20 * π)
prob = ODEProblem(mathieu!, u0, tspan, p)
sol = solve(prob, Tsit5())

Such definition doesn’t produce the correct result, and the linter in VS Code warns about du:

An argument is included in a function signature but not used within its body.Julia(UnusedFunctionArgument)

However, if I define the system like in the tutorial

function mathieu!(du, u, p, t)
du[1] = u[2]
du[2] = # ...
end

everything works fine.

Why does my definition not work? If the u vector had 1000 of components, would I have to write every one of them by hand?

Not just in in-place functions or for vectors. In general, x = … says “attach the label x to ….” If x is an object that already exists, this is different from changing/mutating the object x!

What .= says instead is “do x[i] = … for every i. Unlike x = …, x[i] = ... really does change x!

While du .= A(t, p) * u would work here, use mul!(du, A(t, p), u) to avoid allocation in that operation altogether. You will need to using LinearAlgebra to access that function.

You should be able to write this more simply (and likely with fewer allocations) as return [0 1; -(ω + ε * cos(t)) 0]. Or, if you don’t mind it being a little less “clean” you could write out the operation manually as

function mathieu!(du, u, p, t)
du[1] = u[2]
du[2] = -(p[1] + p[2] * cos(t)) * u[1]
return du
end

Since you’re working with tiny vectors and matrices, you might eventually check out StaticArrays.jl. They can be considerably faster at small sizes. Although it may take a little work to adapt this to it, so only try that once you have this working and if you think a little extra performance is useful.