xMin=0;
xMax=1000;

yMin=0;
yMax=1000;

VxMin=-5000;
VxMax=5000;

VyMin=-5000;
VyMax=5000;

t=100000;

global eps sigma N m radius U

eps=1; %LJ paramater
sigma=100;%LJ parameter
N=20;
m=10;
radius=50;
U=[1/5*xMax;2/5*yMax];

dt=10^-5;

A=zeros(4*N,1);

for i=1:4:4*N-3
    A(i)=xMax*rand();%initial x values
    A(i+1)=yMax*rand();%initial y values
    while norm([A(i);A(i+1)]-U)<=radius
        A(i)=xMax*rand();
        A(i+1)=yMax*rand();
    end    
    A(i+2)=VxMin+(VxMax-VxMin)*rand();%initial Vx values
    A(i+3)=VyMin+(VyMax-VyMin)*rand();%initial Vy values
    
end

for i=1:t
    F1=ode(A);
    F2=ode(A+dt/2*F1);
    F3=ode(A+dt/2*F2);
    F4=ode(A+dt*F3);
    Anew=A+dt/6*(F1+2*F2+2*F3+F4);
    for k=1:4:4*N-3
        %reflective boundary conditions
        if Anew(k)<0 || Anew(k)>xMax
            Anew(k)=A(k);
            Anew(k+2)=-A(k+2);
        end
        if Anew(k+1)<0 || Anew(k+1)>yMax 
            Anew(k+1)=A(k+1);
            Anew(k+3)=-A(k+3);
        end
        %condensation
        if norm([Anew(k);Anew(k+1)]-U)<=radius
            Anew(k+2)=0;
            Anew(k+3)=0;
        end
    end
    A=Anew;
    B=zeros(N,1);
    C=zeros(N,1);
    for j=1:N
        B(j)=A(1+4*(j-1));
        C(j)=A(2+4*(j-1));
    end
    X=U(1)-radius:1:U(1)+radius;
    Y1=U(2)+(radius^2-(X-U(1)).^2).^(1/2);
    Y2=U(2)-(radius^2-(X-U(1)).^2).^(1/2);
    if mod(i+99,100)==0
    plot(X,Y1,'k-')
    axis([xMin xMax yMin yMax])
    hold on
    plot(X,Y2,'k-')
    plot(B,C,'ro')
    pause(0.01)
    hold off
    end
end
