test2.m

tm=10;
dt=0.01;
ta=0:dt:tm;
Lta=length(ta);

d0=0.01; % distance of stability
ss=0.03; % square size
m=1e-3;
k=1; % spring coefficient
fcr=2*sqrt(k/m); % critical friction coeffcient
f=2*fcr; % friction coeficent


r=zeros(2,0);
g=[0; -10]; % gravitation

R=0.1;
R2=R^2;
ph=-0.1; % plate level
% make circle:
for yt=-R:d0:R
    for xt=-R:d0:R
        if (xt^2+yt^2<=R2)
            r=[r [xt;yt]];
        end
    end
end
hr=plot(r(1,:),r(2,:),'k.');
hold on;
hpl=plot([-0.15 0.15],ph*[1 1],'r--');
axis equal;

v=zeros(size(r)); % velocity
N=size(r,2);
Na=1:N;

for tc=2:Lta
    t=ta(tc);
    D=NaN(N,N); % distance matrix
    for n1=1:N-1
        rt=r(:,n1);
        ind2=n1+1:N;
        ft=find((abs(rt(1)-r(1,ind2))<=ss)&(abs(rt(2)-r(2,ind2))<=ss)); % rest points in square
        d=sqrt((rt(1)-r(1,ft)).^2+(rt(2)-r(2,ft)).^2);
        ind2ft=ind2(ft);
        D(n1,ind2ft)=d;
        D(ind2ft,n1)=d;
        
    end
    
    A=(D-d0)/m; % acelerations from edges
    v=v+dt*;
end