Well I'll just leave this here: I'm trying to implement in DiffEq the update rule for a spiking neural network. For now I've implemented with plain Julia and works decently, but I'd like to speed it up a little:

The function simulate just initialize the variables and then runs step! for 10000 times.

The function step! takes as input:
V: a matrix(neurons x time), where values are potentials of individual neurons at a certain time
sd: a vector(neurons), where values represent how long ago a given neuron spiked
w: a matrix(neurons x neurons) of weights (connectivity matrix)
np, ep are parameters
tottime is the total time of simulation

function simulate(tottime, np, sp, ep)
    w = genweights(np, sp)
        V = zeros(np.N, tottime)
        sd = zeros(np.N) .+ 0.2

    V[:, 1:4] .= np.Vᵣ
    V[126:132,1:4] .= 0.

    for t in 3:tottime-1
        step!(V, sd, t, w, np, ep, tottime)
    end

    V[:, 1:tottime], sd
end


function step!(V, sd, t, w, np, ep, tottime)

    V[:, t+1] = dv(view(V, :, t), sd, w, np, ep) # This is the update
    V[view(V,:, t) .> np.Vₜ, t+1] .= 0.0  # This makes the neurons over the threshold spike
    V[view(V,:, t) .== 0.0, t+1] .= np.Vᵣ  # This reset the neurons that spiked to a low potential 

    sd .+= ep.dt  # Here we adjust the spike delays for the next timestep
    sd[view(V,:, t-2) .== 0.0] .= 3*ep.dt  # This sets the spike delays to 3 timesteps for the neurons who spiked 2 timesteps ago, emulating transmission delays or something
    end

end

function dv(v, sd, w, np, ep)
    k = sd .* exp.(-sd ./ np.τₛ) # This is just to save computations
    v.+ ((Eₗₜ .- v ./ np.τₘ ).- (v.-np.Eₑ) .*( w[2] * k) .- (v.-np.Eᵢ) .*( w[1] * k)) .* ep.dt  # This is the differential equation, evaluated with Euler
end

The function step! runs 10000 times, simulating 10 seconds of activity, and it takes ~350ms
I'd like to implement this in DiffEq, but I don't understand how to write a function in the way DiffEq wants it and having it computing the spikes.

Specifically I'm confused by:

• How to keep track of the timestep inside the function
• How I should define the matrix
• What will solve do to my matrix: should I multiply by dt, or will solve do it?

I came up with this, I don't know if it makes any sense. I think there is a main problem in updating sd

# This function initializes all the needed arrays
function DEsimulate(tottime, np, sp, ep)
    alg = FunctionMap{true}()
    w = genweights(np, sp)
    v = zeros(np.N) .+ np.Vᵣ
    sd = zeros(np.N) .+ 0.2

    v[1:5] .= 0.

    p = (w, np, ep)
    prob = DiscreteProblem(update!, [v, sd], (ep.dt, tottime*ep.dt), p)
    sim = solve(prob, alg, dt=ep.dt)
    return sim

end

function update!(dv, u, p, t)
    w, np, ep = p

    v, sd = u
    v[v .== 0.0] .= np.Vᵣ
    v[v .> np.Vₜ] .= 0.0

    k = sd .* exp.(-sd ./ np.τₛ)
    dv[1] += ((Eₗₜ .- v ./ np.τₘ ).- (v.-np.Eₑ) .*( w[2] * k) .- (v.-np.Eᵢ) .*( w[1] * k))
    dv[2][v .== 0.0] .= 2*ep.dt
    dv[2] .+= 1
end

Any kind of hint, also on style or everything else is very well appreciated!