Размеры матриц одной переменной не соответствуют

Занимаюсь решением следующей физической задачи: необходимо рассчитать динамику пучка электронов, двигающихся в поле квадруполя. Для генерации точек использую распределение по Гауссу, и формирую сначала 3 массива данных по координатам x,y,z, а далее формирую массив радиус векторов, каждый элемент которого это массив переменных x,y,z отдельного электрона (такую модель посоветовал руководитель работы). Далее, не вдаваясь в подробности, расчет идет по следующему пути: считается сила Лоренца как векторное произведение, на ее основе определяется изменение импульса, и, соответственно, сам импульс как предыдущее его значение+приращение. И вот тут возникает загвоздка: переменная имеет тип 3х10 double, как и все матричные значения в расчетах, но уже на этапе расчета импульса возникает ошибка несоответствия размеров

Unable to perform assignment because the size of the left side is 1-by-10 and the size of the right side is 3-by-10.

Пробовал менять размеры матриц в разных вариациях, расширять матрицу внутри цикла перед расчетом, однако результата не получил. Ниже прикреплю последний вариант кода, я его регулярно редактирую, и пока дальше этой ошибки не разбирался, возможно подобная проблема возникает не только на этом этапе. Заранее благодарю всех, кто сможет чем-то помочь.

clear
%Generation points
rng('default') % For reproducibility
X=[];
Y=[];
Z=[];
mu = 1;
sigma = 0.01;

for i=1:10
x = random('Normal',mu,sigma);
y=random('Normal',mu,sigma);
z=random('Normal',mu,sigma);
X(end+1)=x;
Y(end+1)=y;
Z(end+1)=z;
R=[X;Y;Z];
end
%Variables
N=10;
u=10;
s=10;
dt=1*10^-6;
q=1.6*10^-19;
eps=8.85*10^-12;
g=15;
gz=zeros(1,u);
Vz=repmat(0.5*10^5,s,u);
z=cumsum(Vz*dt);
V=zeros(3,s);
m=9.1*10^-31;
p=zeros(3,s);
%Calculation
for i=1:N-1
    B=[g*X(i,:); g*Y(i,:); gz];
    F=q*cross(B,V);
    dp=F*dt;
    p(i+1,:)=p(i,:)+dp;
    V(i+1,:)=p(i+1,:)/m;
    dr=V(i+1,:).*dt;
    R(i+1,:)=R(i,:)+dr;
end
R