Turn system of ODE into vectorized form

Hi, I am trying to solve an ode system of ~5000 trajectories each with 3 components {x,y,z} for a total of ~15000 variables. I have tried to vectorize my code to speed it up and it works fine but I have yet to vectorize the ODE itself. Currently I am using y(1) -> y(5000) as the X components y(5001) -> y(10000) as the Y components and y(10001) -> y(15000) as the Z components. How would I go about finalizing this vectorization and will I see an increase in speed? Also, I'm having trouble with running out of memory, will this help with that issue?

function dSdt = eomVect(t,S)
%eom generate equations of motions (VECTORIZED)
%
%xlam0/mu0/delta0 = .00148868809/5.179275814074202/0.005078629803839
%
%xlamf/muf/deltaf = 0.001646784605151/5.127046260784226/.108   0,80
%xlamf/muf/deltaf = 0.001877329529838/2.722334930113119/.961   0,10

nc = .825;
wD = 50*2*pi*nc;
kD = sqrt(2*wD);
dk = .404953;
xlamf = 0.001877329529838;

[kx ky] = meshgrid(-kD:dk:kD,0:dk:kD);
kvec = sqrt(ky.^2+ kx.^2)<=kD;
kxvec = kx(kvec);
kyvec = ky(kvec);

Nspins = size(kxvec,1);

dSdt = zeros(Nspins*3,1);

Sx = S(1:Nspins);
Sy = S(Nspins+1:2*Nspins);
Sz = S(2*Nspins+1:3*Nspins);

%calculate delta_x and delta_y
epsk = kxvec.^2+kyvec.^2;
sqepsk = sqrt(epsk);
DX = sum(-2*dk^2/pi*xlamf*sqepsk.*Sx);
DY = sum(-2*dk^2/pi*xlamf*sqepsk.*Sy);

Bx = -2*sqepsk*DX;
By = -2*sqepsk*DY;
Bz = -epsk;

dSdt(1:Nspins)= Sy.*Bz-Sz.*By; %x
dSdt(Nspins+1:2*Nspins)= Sz.*Bx-Sx.*Bz; %y
dSdt(2*Nspins+1:3*Nspins)= Sx.*By-Sy.*Bx; %z

end