Метод молекулярной динамики и схема Верле

clear;
N = 100; %количество частиц
Tk = 60;  %время окончания расчетов
dt = 0.1; % шаг времени
Nh = Tk/dt; %количество шагов
 
L = 1000; %длина ящика
H = 1000; %ширина
 
%задаем начальные координаты и скорости
for i = 1:N
    x(i) = rand*L;
    y(i) = rand*H;
    Vx(i) = rand;
    Vy(i) = rand;
end;
 
%начальные ускорения частиц равны нулю
for i = 1:N
    ax(i) = 0;
    ay(i) = 0;
end;
 
%главный цикл
for t = 1:Nh
    for i = 1:N 
        x(i) = x(i) + Vx(i)*dt +0.5*ax(i)*dt^2;
        y(i) = y(i) + Vy(i)*dt +0.5*ay(i)*dt^2; 
      
        %проверка граничных условий - они должны быть периодическими
      if x(i)<0
          x(i) = x(i) + L;
      end;
      if x(i)>L
          x(i) = x(i) - L;
      end;
      if y(i)<0
          y(i) = y(i) + H;
      end;
      if y(i)>H
          y(i) = y(i) - H;
      end;
      
      %вычисление скорости со старым ускорением
      Vx(i) = Vx(i) + 0.5*ax(i)*dt;
      Vy(i) = Vy(i) + 0.5*ay(i)*dt;
      
      %%%%%%%%%%%%%%%%%
      % вычисление нового ускорения
      for j = i+1:N
 
          %???????????????
 
      end;
      %%%%%%%%%%%%%%%%%
      
      %вычисление скорости с новым ускорением
      Vx(i) = Vx(i) + 0.5*ax(i)*dt;
      Vy(i) = Vy(i) + 0.5*ay(i)*dt;    
      
        
    end;
    plot(x,y,'*')
end