问在Julia中求解初始条件为二维数组的常微分方程组

EN

http://diffeq.sciml.ai/latest/tutorials/ode_example.html

我正在尝试使用Julia (DifferentialEquations.jl)中的常微分方程求解器来求解n个相互作用的粒子系统。假设系统是二维的，每个粒子的运动方程是由它相对于时间的位置的二阶常微分方程来描述的。然后每个粒子需要四个变量，两个用于位置，两个用于速度。然后需要声明4n个变量。有没有一种方法可以推广，这样就不需要逐个列出所有4n个方程？

例如：

http://diffeq.sciml.ai/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1

我尝试将上面链接中的洛伦兹方程修改为n个粒子(这是一个非常非常粗糙的尝试，因为我实际上不知道如何做到这一点)，尝试将u和du扩展到2D阵列。

using DifferentialEquations
using Plots
n = 4
function lorenz(du,u,p,t,i)
    du[i,1] = 10.0*(u[i,2]-u[i,1])*sum(u[1:n,1])
    du[i,2] = (u[i,1]*(28.0-u[i,3]) - u[i,2])*sum(u[1:n,1])
    du[i,3] = (u[i,1]*u[i,2] - (8/3)*u[i,3])*sum(u[1:n,1])
end

u0 = hcat([1.0;0.0;0.0], [0.0;1.0;0.0], [0.0;0.0;1.0])
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob)

这并不令人惊讶，但我希望您能理解我想要做的事情。有没有办法让ODE求解器将u作为2D数组来求解(或者其他可以实现类似目的的方法)？