Не получаются правильные расчеты

Задача следующая:
10 частиц распределены вдоль линии по Z, т.е. по направлению движения
Расстояние между частицами одинаково.
Q = 10 нКл – заряд каждой частицы.
Масса каждой частицы равна массе электрона.
L = 0,004 см – длина последовательности частиц.
N= 10 - кол-во частиц,
W = 10 МеВ – начальная полная энергия частиц.
Частицы начинают двигаться вдоль положительного направления оси Z
Найти распределение энергии между частицами после прохождения 1 см пути, учитывая только Кулоновское взаимодействие.* Для решения задачи я написал следующий код в Matlab (версия R2020b)

clc
clear
N=10; 
dt=1*10^-12; 
q=10^-9;
eps0=8.85*10^-12;
c=3*10^8;
pi=3.1415926;
m=9.1*10^-31;
k=1/(4*pi*eps0);
W=[1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13];
R=[0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0 0; 4*10^-5 8*10^-5 12*10^-5 16*10^-5 20*10^-5 24*10^-5 28*10^-5 32*10^-5 36*10^-5 4*10^-4];
dp=zeros(3,N);
dr=zeros(3,N);
V=zeros(3,N);
F_Q_2=zeros(3,N,N);
F_Q=zeros(3,N);
F_Q2=zeros(3,N);
Z=zeros(1,N);
p=zeros(3,N);
X=[];
Y=[];
Z=[];
R=vpa(R); 
while R(3,10)<4*10^-4+10^-2 
    F_Q = ForceQ(F_Q, R,N,k,q);
  for i=1:N
    p(3,i)=sqrt((W(1,i)^2/c^2)-(m^2*c^2));
    V(3,i)=((c^2/W(1,i))*p(3,i));
    dp=F_Q(1:3,i)*dt;
    p(1:3,i)=p(1:3,i)+dp;
    pn=sqrt(p(1,i).^2+p(2,i).^2+p(3,i).^2);
    W(1,i)=c*sqrt(pn^2+m^2*c^2);
    V(1:3,i)=(c^2/W(1,i))*p(1:3,i);
    dr(1:3,i)=V(1:3,i).*dt;
    R(1:3,i)=R(1:3,i)+dr(1:3,i);
    R = real(R);
    Z(1,i) = R(3,i);
   end
 end
colormap hot
plot(Z(1:N),W(1:N), '.');
ylabel('\itW, Дж','fontsize',14)
xlabel('\itZ, м','fontsize',14)

function F_Q = ForceQ(F_Q2, R,N,k,q)
R = circshift(R,1,2);
for count=1:N
    F_Q_2=zeros(3,N,N);
    R = circshift(R,-1,2);
    for var=2:N
        r=sqrt((R(1,var)-R(1,1))^2+(R(2,var)-R(2,1))^2+(R(3,var)-R(3,1))^2);
        force=(k*q^2)/r^2;
        alpha=asin(real(R(3,var)-R(3,1))/r);
        rxy=r*cos(alpha);
        gamma=asin(real(R(2,var)-R(2,1))/rxy);
        gamma2=asin(real(R(1,var)-R(1,1))/rxy);

        F_Q_2(1,1,var)=F_Q2(1,count)-(force*1*sin(gamma2));
        F_Q_2(2,1,var)=F_Q2(2,count)-(force*1*sin(gamma));
        F_Q_2(3,1,var)=F_Q2(3,count)-(force*sin(alpha)); 
    end
    F_Q1 = sum(F_Q_2,3);
    F_Q(:,count) = F_Q1(:,1);
 end
 end

Опуская физику, чисто математически мне нужно построить зависимость W от Z, причем частицы должны менять координаты таким образом, что первые 5 двигаются в отрицательном направлении, оставшиеся 5 - в положительном. В корректности формул я уверен, однако расчеты получаются не совсем верными. Проверял в debug, первые два прохода в цикле while считаются корректно, на 3 проход уже при вызове функции F_Q на выходе значения численно верные, но знаки у 5 и 6 частиц инверсно меняются (должно быть -число у 5 и +число у 6, а на выходе наоборот). Причем функцию F_Q проверял в отдельном файле с циклом for, считается корректно. Соответственно, из-за этого портятся все остальные расчеты координат и на выходе к концу цикла while получается совершенно не то, что должно быть. Пробовал также использовать другой вариант расчета F_Q без использования функции, через циклы, однако там проблема встает при расчете проекций вектора. Код вот (прикрепляю уже сам цикл, тк начальные условия до цикла те же):

    while R(3,10)<4*10^-2
for i=1:N
    p(3,i)=sqrt((W(1,i)^2/c^2)-(m^2*c^2));
    V(3,i)=((c^2/W(1,i))*p(3,i));
 for var=i+1:N
    r=sqrt((R(1,var)-R(1,i))^2+(R(2,var)-R(2,i))^2+(R(3,var)-R(3,i))^2);
    force=(k*q^2)/r^2;
    alpha=asin(abs(R(3,var)-R(3,i))/r);
    rxy=r*cos(alpha);
    gamma=asin(abs(R(2,var)-R(2,i))/rxy);
    F_Qx(1,i)=F_Qx(1,i)+(force*cos(alpha)*cos(gamma));
    F_Qy(1,i)=F_Qy(1,i)+(force*cos(alpha)*sin(gamma));
    F_Qz(1,i)=F_Qz(1,i)+(force*cos(alpha));
    F_Qx(1,var)=F_Qx(1,var)-(force*cos(alpha)*cos(gamma));
    F_Qy(1,var)=F_Qy(1,var)-(force*cos(alpha)*sin(gamma));
    F_Qz(1,var)=F_Qz(1,var)-(force*cos(alpha));
    F_Q(1:3,i)=[F_Qx(1,i); F_Qy(1,i); F_Qz(1,i)];
    F_Q(1:3,var)=[F_Qx(1,var); F_Qy(1,var); F_Qz(1,var)];
    end
end
for i=1:N
 F=F_Q(1:3,i);
 dp=F*dt;
 p(1:3,i)=p(1:3,i)+dp;
 pn=sqrt(p(1,i).^2+p(2,i).^2+p(3,i).^2);
 W(1,i)=c*sqrt(pn^2+m^2*c^2);
 V(1:3,i)=(c^2/W(1,i))*p(1:3,i);
 dr(1:3,i)=V(1:3,i).*dt;
 R(1:3,i)=R(1:3,i)+dr(1:3,i);
 R = real(R);
 Z = R(1:N);
 end
end
colormap hot
plot(Z(1:N),W(1:N), '.');
ylabel('\itW, Дж','fontsize',14)
xlabel('\itZ, м','fontsize',14)

Я не понимаю в чем причина неправильных расчетов, мои предположения в том, что либо я неправильно оформил порядок в цикле, либо расчет портит circshift (эту функцию я использовал в другой похожей программе, расчет F_Q велся аналогично, и в целом там все работало корректно).