Создать анимацию движения частиц друг относительно друга

clear;
clc;
clf;
 
n = 4;  % кол-во точек
vR = 1; % радиус Ван-дер-Ваальса
 
P = randn(2,n);
x = P(1,:);
y = P(2,:);
dx = 5;
ang=[];
for q = 1:50 % отображение 50 шагов анимации
    
    for i =1:n
        j = 1:n;
        % Расчет результирующих векторов
        dx(i) =  - sum( x(i) - x(j) );
        dy(i) =  - sum( y(i) - y(j) );
        % Углы между цетрами точек
        ang(i) = - atan( dy(i)/dx(i) );
    end
    
    % Корретировка вектороа  учетом области Ван-дер-Ваальса
    dx = 0.1*(dx + vR.*cos(ang))/n;
    dy = 0.1*(dy + vR.*sin(ang))/n;
    % Получение новых координат
    x = x + dx;
    y = y + dy;
    
    % Прорисовка
        plot(x,y,'.r')
        hold on
        axis equal
        grid off
        title (['Frame: ',num2str(q)]);
     
    pause(0.2)
   
end

clear;
N=2;
 
if N==2,
   dt=1e-19;
   X=3;
   Y=3;
elseif N==3,
   dt=1e-21;
   X=5;
   Y=5;
elseif N==50,
   dt=1e-15;
   X=100;
   Y=100;
end;
 
T=100;
m=1e-15;
a=1e+20;
b=1e+20;
 
for ii=1:N,
    x(ii)=X*rand(1);
    y(ii)=Y*rand(1);
end
 
%x=[2.1 1.05];
%y=[0 1.1*cos(pi/6)];
%x=[0 2.1 1.05];
%y=[0 0 1.1*cos(pi/6)];
 
 
plot(x,y,'ro')
pause;
Vx(1:N)=0;
Vy(1:N)=0;
I=0;
for t=0:dt:T,
I=I+1;
    for ii=1:N,
    Fx=0;
    Fy=0;
    for jj=1:N,
        if jj~=ii,
        R=((x(ii)-x(jj)).^2+(y(ii)-y(jj)).^2).^(0.5);
        alf=atan(y(jj)-y(ii))/(x(jj)-x(ii));
        alf=abs(alf);
        F=-a.*R.^(-13)+b.*R.^(-7);
        if (x(jj)-x(ii))>=0 && (y(jj)-y(ii))>=0,
            Fx=Fx+F*cos(alf);
            Fy=Fy+F*sin(alf);
        elseif (x(jj)-x(ii))<0 && (y(jj)-y(ii))>=0,  
            Fx=Fx-F*cos(alf);
            Fy=Fy+F*sin(alf);
        elseif (x(jj)-x(ii))<0 && (y(jj)-y(ii))<0,
            Fx=Fx-F*cos(alf);
            Fy=Fy-F*sin(alf);
        elseif (x(jj)-x(ii))>=0 && (y(jj)-y(ii))<0,
            Fx=Fx+F*cos(alf);
            Fy=Fy-F*sin(alf);
        end
        end
    end
    Ax=Fx/m;
    Ay=Fy/m;
    Vx(ii)=Vx(ii)+Ax*dt;
    Vy(ii)=Vy(ii)+Ay*dt;
    dx(ii)=Vx(ii)*dt;
    dy(ii)=Vy(ii)*dt;
end
x=x+dx;
y=y+dy;
r(I)=((x(1)-x(2))^2+(y(1)-y(2))^2)^0.5;
v(I)=((Vx(1)-Vx(2))^2+(Vy(1)-Vy(2))^2)^0.5;
TT(I)=t;
pause(1e-10)
plot(x,y,'ro')
 
XX(I,:)=x;
YY(I,:)=y;
end