How to make the balls in this program organize and settle at the top of the box instead of leaving it

function y = atoms(x)

clc
n = 25;  % no. of atoms
dt = 0.1;   %  time
ma=1;
mb=1;
x0=[0 0 0 0];
l=200;
p = rand(n,2)*0.5*l;
v = rand(n,2)+.5;
r=0.5;
for tt=1:5000

      plot(p(:,1),p(:,2),'O', 'MarkerFaceColor', 'b');
      axis([0  l  0  l])
      drawnow;

      for i = 1 : n
          if p(i, 1) < 0
              v(i, 1) =- v(i, 1);
          end
          if p(i, 1) > l
              v(i, 1) = -v(i, 1);
          end
          if p(i, 2) < 0
              v(i, 2) = v(i, 2);
          end

      end

      for i = 1 : n
          for j = i + 1 : n
              d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);
              if d < r
                  va = sqrt(v(i,1)^2 + v(i,2)^2);
                  vb = sqrt(v(j,1)^2 + v(j,2)^2);
                  gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1)));
                  alpha = atan2(v(i,1),v(i,2)) - gama;
                  beta  = atan2(v(j,1),v(j,2)) - gama;
                  [x,fval ] = fsolve(@gh,x0);

                  v(i,1)=x(1)*cos(x(3)+gama);
                  v(i,2)=x(1)*sin(x(3)+gama);
                  v(j,1)=x(2)*cos(x(4)+gama);
                  v(j,2)=x(2)*sin(x(4)+gama);

      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;

              end
          end
      end

      p=p+v*dt;

end

      function y = gh(x)
          e = 1;
          y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
               va*sin(alpha) - x(1)*sin(x(3));
               vb*sin(beta)  - x(2)*sin(x(4));
               x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];

      end

end

Answer 1

[KSSV]

    function y = atoms()    
    n = 25;  % no. of atoms
    dt = 0.1;   %  time
    ma=1;
    mb=1;
    x0=[0 0 0 0];
    l=200;
    p = rand(n,2)*0.5*l;
    v = rand(n,2)+.5;
    r=0.5;

    for tt=1:5000
        idx = abs(p(:,2)-l)<=1 ;  % added this line         
        p(idx,2) = l ;            % added this line 
        plot(p(:,1),p(:,2),'O', 'MarkerFaceColor', 'b');
        axis([0  l  0  l])
        drawnow;
        for i = 1 : n
            if p(i, 1) < 0
                v(i, 1) =- v(i, 1);
            end
            if p(i, 1) > l
                v(i, 1) = -v(i, 1);
            end
            if p(i, 2) < 0
                v(i, 2) = v(i, 2);
            end
        end
        for i = 1 : n
            for j = i + 1 : n
                d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);
                if d < r
                    va = sqrt(v(i,1)^2 + v(i,2)^2);
                    vb = sqrt(v(j,1)^2 + v(j,2)^2);
                    gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1)));
                    alpha = atan2(v(i,1),v(i,2)) - gama;
                    beta  = atan2(v(j,1),v(j,2)) - gama;
                    [x,fval ] = fsolve(@gh,x0);
                    v(i,1)=x(1)*cos(x(3)+gama);
                    v(i,2)=x(1)*sin(x(3)+gama);
                    v(j,1)=x(2)*cos(x(4)+gama);
                    v(j,2)=x(2)*sin(x(4)+gama);
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                end
            end
        end
        p=p+v*dt;

    end
        function y = gh(x)
            e = 1;
            y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
                va*sin(alpha) - x(1)*sin(x(3));
                vb*sin(beta)  - x(2)*sin(x(4));
                x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];
        end
    end

[zina shadidi]

Thank you very much KSSV . I really appreciate this

[KSSV]

Appreciating is accepting and voting the answer :)

Answer 2

[zina shadidi]

Hi KSSv

Let me ask you please .what if we want to line the balls in more than one row

[KSSV]

Can you elaborate?

Answer 3

[zina shadidi]

The two lines that you added makes The balls lined up at the top of the figure in a row

Now if we want the balls to line up in more than one row, ???

[zina shadidi]

the above image shows the balls lined in one image