m1 = 1;
m2 = 1;
k = 30;
l1 = 1;
l2 = 3;
l3 = 1;
d = 0.3;

x1 = 2
v1 = 0;
x2 = 3;
v2 = 0;

T1 = 0:0.005:10;
T2 = 0:0.008:10;
delay = 0.2;

X1 = [x1];
X2 = [x2];
V1 = [v1];
V2 = [v2];

i1=2;
i2=2;

while i1<=length(T1) | i2<=length(T2)
    x1 = X1(end);
    v1 = V1(end);
    x2 = X2(end);
    v2 = V2(end);
    if i2>length(T2) | (i1<=length(T1) & T1(i1)+delay < T2(i2))
        dt=T1(i1)-T1(i1-1);
        f11 = -k * (x1-l1)/l1 - d*v1;
        f12 = k * ((x2 - x1)-l2)/l2 + d*(v2-v1);
        a1 = (f11+f12)/m1;
        %x1 = x1 + dt*v1 + 0.5*dt*dt*a1;
        x1 = x1 + dt*v1;
        v1 = v1 + dt*a1;
        X1 = [X1 x1];
        V1 = [V1 v1];
        i1 = i1+1;
    else
        dt=T2(i2)-T2(i2-1);
        f21 = -k * ((x2 - x1)-l2)/l2 - d*(v2-v1);
        f22 = k * (l1+l2+l3-x2-l3)/l3 - d*v2;
        a2 = (f21+f22)/m2;
        %x2 = x2 + dt*v2 + 0.5*dt*dt*a2;
        x2 = x2 + dt*v2;
        v2 = v2 + dt*a2;
        X2 = [X2 x2];
        V2 = [V2 v2];
        i2 = i2+1;
    end
end

[Ta, X1a, V1a, X2a, V2a] = sync(0.001);
[Ts, X1s, V1s, X2s, V2s] = sync(0.005);
close all
figure
plot(T1,X1,'g');
hold on
plot(T2,X2,'b');
plot(Ts,X1s,'r');
plot(Ts,X2s,'r');
%plot(Ta,X1a,'g');
%plot(Ta,X2a,'g');
