Hello,

I normally deal with dynamics on networks, so I have giant sparse systems of coupled ODE’s. The problem I have at the moment is that for a simple discrete map, a naive implementation is significantly faster than solving the DiscreteProblem provided by DIfferentialEquations.jl. For example, the code below is a MWE of what my problems look like,

```
using DifferentialEquations
using Random
using SparseArrays
# parametes
rng = MersenneTwister(0)
N = 80000
M = 150000
T = 10
β = 0.4
A = sprand(rng, N, M, 0.0001, (r, x) -> rand(r, [1], x))
p = (β, A)
###
function map!(unp1, un, p, t)
β, A = p
unp1 .= β .* un .* A
end
function naive_test(u0, tspan, p)
unp1 = similar(u0)
for i in tspan[1]:tspan[2]
map!(unp1, u0, p, i)
end
return unp1
end
# test
u0 = sprand(N, M, 0.0001)
prob = DiscreteProblem(map!, u0, (0, T), p)
@time solve(prob; save_on=false);
@time naive_test(u0, (0, T), p);
```

For the `solve`

function I turn off saving of intermediate results to better match the naive test and I get the following,

```
julia> @time solve(prob; save_on=false);
205.226954 seconds (16.25 M allocations: 1.241 GiB, 0.10% gc time, 1.79% compilation time)
```

and for the `naive_test`

implementation I get

```
julia> @time naive_test(u0, (0, T), p);
0.142919 seconds (16.55 k allocations: 20.360 MiB, 7.48% compilation time)
```

I was wondering if I’m missing something and this sort of performance is expected, I know there is a significant amount of machinery behind diffeq and it might not suit sparse arrays. Any help or explanation would be appreciated.

Using DifferentialEquations v7.1.0 and julia 1.7.1.

Rory.