julia differential equations for system of many equations isn't working

I'm changing my programming language to julia due to it's better performance for numerical calculations and differential equations compared to python (I used mostly scipy solve_ivp and odeint but they turned out being too slow for my problems).

For testing purposes I did a simple translation from my old code in python to a new code in julia and I'm pretty sure both codes should have the same numerical results. The problem seems to appear when I try to use the ODE solver from julia and I can't figure out why (maybe because there is too many equations?)

My physics problem is a simulation of a beam of charged planes. For the integration over the time, I use a 2n-dimensional vector, with the [1:n] coordinates corresponding to the positions of the planes and the [n+1:2n] positions correspondind to the particles speeds. So the derivative function is 2n-dimensional as well and the [1:n] time derivatives corresponds to the [n+1:2n] coordinates of the initial vector and the [n+1:2n] derivatives corresponds to a positional term and an index term (see codes below).

Starting by my working code in python:

def deriv(t,y):

        p= len(y)
        dydt = np.zeros(int(p))

        dydt[0:int(p/2)] = y[int(p/2):int(p)] #derivative of the positions

        k = index(y[0:int(p/2)])
        signal=abs(y[0:int(p/2)])/y[0:int(p/2)]

        dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*signal 
        #derivative of the speeds

        deriv = dydt[0:int(p)]

        return deriv

With the following integration

sol = solve_ivp(deriv, (time,time+ delta_t), initial_phase, rtol = 0.00001)

and the resulting graphic of the phase-space looks like this:

Few time expected evolution:

Then we have my code in Julia:

function dydt(du, u,p,t)
        n = floor(Int,length(u)/2)

        k = sortperm(abs.(u[1:n]))
        k = k[k] #k ends up beeing the index equivalent

        du[1:n] = u[n+1: 2n] #derivative of the positions

        du[n+1: 2n] = (1/2n)*(-ones(n) + 2*k).*(abs.(u[1:n])./u[1:n])- (u[1:n]) 
        #derivative of the speeds
 end

With the following integration:

using DifferentialEquations
p = 0  #there is no p parameters in my code but they are defined in the tutorials anyway
H = 1.5*rand(2000)
L = (2/3)*(rand(2000) - 0.5*ones(2000))
D = append!(H,L) #Initial state, similar to the python example 
tspan = (0.0,1.0)
prob = ODEProblem(dydt,D,tspan,p,reltol=1e-4,abstol=1e-4)
sol = solve(prob)

And the result is assimetrical as we can see below:

I hope anybody can help. I'll answer any asked major question about the physics envolved if necessary.