forwarddiff with in-place functions

fd.jl

using ForwardDiff

const T = Float64 # the datatype used, probably Float64
const dof = 4

fn!(;T_::Type=T, dof_=dof) = zeros(T_,dof_)
function fn!(res, y)#, dydt)
    @show typeof(res), typeof(y)

    for i=1:dof
        res[i] = 3*y[i]#+dydt[i]
    end
    nothing
end

const res = fn!()
function fn_ex(y)
    fn!(res, y)
    return res
end


ad_jacobian(F) = (t,y)->ForwardDiff.jacobian(F,y)
const out = fn!()
ad_jacobian2(F!) = (t,y)->ForwardDiff.jacobian(F!,out,y)
ad_jacobian2!(F!) = (outJ,t,y)->ForwardDiff.jacobian!(outJ,F!,out,y)

J = ad_jacobian(fn_ex)
J! = ad_jacobian2(fn!)
J!! = ad_jacobian2!(fn!)

ic = zeros(T,dof)


typeof(J(0.0,ic))
typeof(J!(0.0,ic))

const outJ = similar(out, length(out), length(out))
J!!(outJ,0,ic)
typeof(outJ)

You are passing a function with side effects to jacobian, while it probably expects that each evaluation is independent.