correct definition of a differential equation for julia/DifferentialEquations

Question

I have two variants of a code to solve an ODE system y'=A*y+b. This one doesn't work (it runs but gives the wrong result, namely it returns the initial data u0 as a constant solution):

using LinearAlgebra
using DifferentialEquations

function RHSfun(du,u,p,t)
    du= p[1]*u + p[2]
    return du
end


A=[-1 0 0; 0 -1 0; 0 0 -1]

b=[1 ; 2 ; 3 ]

parm = [A,b]


 tspan=(0.0,5.0)
 u0=[0.0; 0.0 ; 0.0]
 prob=ODEProblem(RHSfun,u0,tspan,parm)
 sol=solve(prob)

but if I delete the argument du from RHSfun it works, which seems counter to the documentation (and I have solved ODEs defined with the du argument before).

I am not asking for other ways of writing the problem (i know of several working versions), but want to know why the above does not work. I have solved ODEs with the du argument before, as this is what most tutorial say to do.

Answer 1

If you use the in-place form, you have to make sure you're modifying the the du array, like:

function RHSfun(du,u,p,t)
    du .= p[1]*u + p[2]
    return nothing
end

or more efficiently you can use mul!. Notice that what you return doesn't matter. If you want to use the allocating form, then do not specify a mutating argument:

function RHSfun(u,p,t)
    du = p[1]*u + p[2]
    return du
end