My code solves the Monte Carlo problem of a large ensemble of particles interacting with field. I use the Newton’s Equations to evolve the system. Numerically, ordinary equation is needed. Some pieces of my code are as follows:

```
function Eq1(t, y)
(x, y, px, py) = y
r = sqrt( x^2 + y^2 + s1 )
x_prime = px
y_prime = py
px_prime = -2x/r^3
py_prime = -2y/r^3
[x_prime; y_prime; px_prime; py_prime]
end
function Init2()
x = Array(Float64, ne)
y = similar(x)
px = similar(x)
py = similar(x)
rand1 = rand(ne)
x .= R_init * sin(2pi.*rand1)
y .= R_init * cos(2pi.*rand1)
rand1 = rand(ne)
px .= p_init * sin(2pi.*rand1)
py .= p_init * cos(2pi.*rand1)
println(R_init)
println(p_init)
println(s1)
for i in 1:ne
tout, yout = ode45(Eq1, [x[i], y[i], px[i], py[i]], [0.0; 100.0])#; points =:specified)
x[i] = yout[2][1]
y[i] = yout[2][2]
px[i] = yout[2][3]
py[i] = yout[2][4]
if mod(i, 1000) == 0
println("i=$(i)")
end
end
```

Whe I set the parameters as

```
R0_init=0.5
p0_init=0.7869
s1 = 0.5
ne = 5000
```

the @time Inite() command give a result as:

i=1000

i=2000

i=3000

i=4000

i=5000

14.935459 seconds (267.48 M allocations: 8.997 GB, 4.63% gc time)

I have solved similar problem with Fortran and Matlab, both of them are much more efficient. I guess I have wrongly used the package, but I can only find very little details in the package github page. Can anyone tell me what the problem is and how can i fix it?