I am trying to implement a system of two ODEs. The following works:

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
du[1] = -αc
du[2] = -βce
du = [dc, dce]
end

u0 = [3.0, 7.0]
tspan = (0., 2.)
p = [1.0, 1.]

prob = ODEProblem(diffusion!, u0, tspan, p)
sol = solve(prob, Tsit5())

But the does not: notice the difference in the diffusion! function.

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
dc = -αc
dce = -βce
du = [dc, dce]
end

The second version leads to a solution always equal to the initial conditions, while the first case leads to the proper solution. I do not understand why the line

Instead updating the contents of du, with the line du = [dc, dce] you are creating a new vector and setting du to point (binding) to that new vector. The underlying DiffEq code is still using the old one and so does not see any change.

An alternative syntax that would work is broadcasting du .= [dc, dce] since this effectively expands out into the syntax in your first example where you mutate the underlying vector.

Is here a typo here [first example]… just before end, you define du = [dc,dce], but dc and dce have not been defined…I assume that your example should be:

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
du[1] = -α c
du[2] = -β ce
#du = [dc, dce]
end

Thank you!! That is what I was missing. Coming from the world of C++ and Python, it will take me a while to get used to this approach using broadcasting.