I am trying to numerically integrate a set of master equations with DifferentialEquations.jl

I have defined the differential equation like:

```
function voter_model_ame!(dG,G,t,p)
N,U = size(G)
for n_prime=1:N,u_prime=1:U
u = u_prime-1
n = n_prime-1
dG[n_prime,u_prime] = 0
if n == 0
continue
end
if u > n
continue
end
(u < n) && (dG[n_prime,u_prime] -= G[n_prime,u_prime]*(n-u)*(u/n))#upflip outflux
(u > 0) && (dG[n_prime,u_prime] -= G[n_prime,u_prime]*(u)*((n-u)/n))#downflip outflux
(u < n) && (dG[n_prime,u_prime]+= G[n_prime,u_prime+1]*((u+1)*((n-u-1)/n)))#inwards downflip flux
(u > 0) && (dG[n_prime,u_prime]+= G[n_prime,u_prime-1]*((n-u+1)*((u-1)/n)))#inwards upflip flux
end
return dG
end
```

This is a system of N^2 differential equations.

I can set the parameter N = 500 and I try to numerically integrate like this:

```
t_max = 10
p = [0]
u₀ = SAME.initialize_u0(N_individuals,p_dist)#creates initial conditions, Matrix{BigFloat} of size NxN
t_span = (0.,t_max)
prob = ODEProblem(voter_model_ame!,u₀,t_span,p)
sol = solve(prob,Tsit5(),saveat=1)
```

The code take ~3 seconds for N = 50 and upwards of 10min to run with N = 500.

Here are the @time return for N = 50

`4.603619 seconds (77.25 M allocations: 3.787 GiB, 8.55% gc time, 46.87% compilation time)`

Now this seems to be very slow and takes up an insane amount of memory.

I used the in-place style function instead of allocating a new array every time step as per the optimization instructions in the DifferentialEquations.jl I can’t use static arrays because the size of my arrays is too large.

I have a feeling I am doing something horrendously wrong but I don’t know what. There must be an efficient way to solve large systems of differential equations. What am I missing?