When setting up the test for the hemodynamic response model I noticed that I wasn't able to reproduce the numerical solution of the respective differential equations with the Julia implementation of the following equation that spm_int_J allegedly computes
dx = (expm(dt*J) - I)*inv(J)*fx
This was done with ExponentialUtilities.jl
Δx = (exponential!(dt*J[:,:,2]) - I)*inv(J[:,:,2])*reshape(dx[:,:,2],nd-nr)
Putting this on the back burner for now but dropping this Issue here to not forget to check this more thoroughly later.