bacterial.jl

N = 10
λ=-0.1; tEnd=1000.
tThreshold = 0.1
u0 = zeros(4 * N)
for i in 1 : N
     u0[2*i - 1] = -5. + 10.*rand()
     u0[2*i]     = -5. + 10.*rand()
     u0[2*N + 2*i - 1] = 0.
     u0[2*N + 2*i]     = 1.
end
f0= function (t, u, p, du)
    N, λ = Int(p[1]), p[2]
    @inbounds for k = 1 : N
        du[2*N + 2*k - 1] = -u[2*N + 2*k]
        du[2*N + 2*k]     =  u[2*N + 2*k - 1]
        fx, fy = 0.0, 0.0
        for i = 1 : N
            if i != k
                rx = u[2*k-1] - u[2*i-1]
                ry = u[2*k] - u[2*i]
                nx = u[2*N+2*i-1]
                ny = u[2*N+2*i]
                a = rx*nx + ry*ny
                rinv = 1. / √(rx*rx + ry*ry)
                rinv3 = rinv * rinv * rinv
                rinv5 = rinv3 * rinv * rinv
                rinv13 = rinv5 * rinv5 * rinv3
                fx -= rx * rinv3
                fy -= ry * rinv3
                fx += 3*a^2 * rx * rinv5
                fy += 3*a^2 * ry * rinv5
                fx -= rx * rinv13
                fy -= ry * rinv13
            end
            du[2*k - 1] = λ * fx + u[2*N + 2*k - 1]
            du[2*k] =     λ * fy + u[2*N + 2*k]
        end
    end
end
using DifferentialEquations, ODEInterfaceDiffEq
f = ParameterizedFunction(f0, [N, λ])
eq = ODEProblem(f, u0, (0., tEnd))
sol = solve(eq,Rosenbrock23())
sol = solve(eq,CVODE_BDF())
sol = solve(eq,radau())
using Plots; plot(sol,vars=[1],lw=1)