Here is the code:
%%Elastic Collision
clc 
clearvars
close all
%-----------------------
N=5;                  % number of particles
L=10;                   %length of the box
dt=0.0003; 
%X=zeros(2,N); 
%V=zeros(2,N); 
%R=0*rand(1,N)+0.2;   %radius of each particle
%M=2*R; % mass of each particle
theta=0:2*pi/20:2*pi;
nstep=20000;               % number of time steps or duration
pit=200;                   % number of iterations between two plots or speed
%%---------- manually initilizing position adn velocity
 %position
 X=[1 5 0 3 7;...     % x pos    %need to keep 2 row and only expand number of columns
      3 2 4 5 -8];        % y pos
 %velocity
 V=[3 5 2 -10 -3;...       % x
      0 -5 -2 -3 -2];      % y
 %Radius of each particule
 R=[0.6 0.6 0.6 0.6 0.6] ;
 %Mass of each particle
 M=[1 1 1 1 1];
 %
par1=V(:,1);
par2=V(:,2);
par3=V(:,3);
%auxiliary variables  
  Ax=zeros(N,N);
  Ay=zeros(N,N);
  Ar=zeros(N,N);
  for i=1:N
    for j=i:N
        Ar(i,j)=R(1,i)+R(1,j);
    end
  end
 hold on
 axis equal
 for i=1:N       
  x=X(1,i)+R(1,i)*cos(theta);
  y=X(2,i)+R(1,i)*sin(theta);
  plot(x,y,'color','b')
 end
 plotdomain(L)
 grid on
 grid minor
 drawnow
 pause(1)
for T=0:nstep
    X=X+dt*V;
    %check if particules collided with eachother
    for i=1:N
        Ax(i,:)=X(1,i)-X(1,:);
        Ay(i,:)=X(2,i)-X(2,:);
    end
    Ax=triu(Ax);
    Ay=triu(Ay);
    Nrm=(Ax.^2+Ay.^2).^(0.5)-Ar-10^-5; % calculate the distance between each two particles
    Nrm=triu(Nrm(1:N-1,2:N));
    [row,col]=find(Nrm<0); % find the particles that collided 
    l=length(row);
    %if particles collided calculate the new velocities
    if l~=0
        col=col+1;  
        for t=1:l
            i=row(t,1);
            j=col(t,1);
            C1=(2*M(1,j)/(M(1,i)+M(1,j)))*dot(V(:,i)-V(:,j),X(:,i)-X(:,j))/(norm(X(:,i)-X(:,j))^2);
            C2=(2*M(1,i)/(M(1,i)+M(1,j)))*dot(V(:,j)-V(:,i),X(:,j)-X(:,i))/(norm(X(:,i)-X(:,j))^2);
            V(:,i)=V(:,i)-C1*(X(:,i)-X(:,j));
            V(:,j)=V(:,j)-C2*(X(:,j)-X(:,i));
            % change color upon collision % add here??
        end
    end
     for i=1:N
        %check if particules collided with side walls
        if X(1,i)+R(1,i)>=L-10^-5 || X(1,i)-R(1,i)<=-L+10^-5
            V(1,i)=-V(1,i);
        end
        %%check if particules collided with top and bottom walls
        if X(2,i)+R(1,i)>=L-10^-5 || X(2,i)-R(1,i)<=-L+10^-5
            V(2,i)=-V(2,i);
        end 
     end  
     %plot
     if T==pit*ceil(T/pit)
         clf
         hold on
         axis equal
         for i=1:N      
              x=X(1,i)+R(1,i)*cos(theta);
              y=X(2,i)+R(1,i)*sin(theta);
              plot(x,y,'color','b')
         end
        plotdomain(L)
        title(T*dt)
        grid on
        grid minor
        drawnow
        pause(0.01)
     end
end
function plotdomain(L)
    plot([-L,-L],[-L,L],'color',[0 0 0])
    plot([L,L],[-L,L],'color',[0 0 0])
    plot([-L,L],[-L,-L],'color',[0 0 0])
    plot([-L,L],[L,L],'color',[0 0 0])
    plot(-L-0.1,-L-0.1,'color',[1 1 1])
    plot(L+0.1,L+0.1,'color',[1 1 1])
end
