function anime_os(t, v, m1, m2, L1, L2)
  v1=v(1,:);
  v2=v(2,:);

  xr=[ 0 .5 1.5 2.5 3.5 4 4 4];
  xr=xr/4;

  yr=[ 0 1 -1 1 -1 0 1 -1];

  ym=[1 -1 -1 1 1];
  xm=[0 0 1 1 0];

  lm1=2;
  lm2=2*m2/m1;
  clf
  drawlater
      plot(xr*(L1+v1(1)),yr,'linewidth',3);
      r1=get(gce(),'children')
      plot((L1+v1(1))+lm1+xr*(L2+v1(1)-v2(1)),yr,'linewidth',3);
      r2=get(gce(),'children')
      plot(xm*lm1+L1+v1(1),ym,'linewidth',3);
      ma1=get(gce(),'children')
      plot(xm*lm2+L1+L2+v2(1)+lm1,ym,'linewidth',3);
      ma2=get(gce(),'children')
      set(gca(),'isoview','on','auto_scale','off','box','off');
      a=get("current_axes");//get the handle 
      a.data_bounds=[0,-1;22,1];
  drawnow

  for i=1:length(t)
      sleep(50); // <==
	  drawlater
	  set(r1,'data',[xr*(L1+v1(i));yr]');
	  set(r2,'data',[L1+v1(i)+lm1+xr*(L2-v1(i)+v2(i));yr]');
	  set(ma1,'data',[xm*lm1+L1+v1(i);ym]');
	  set(ma2,'data',[xm*lm2+L1+L2+v2(i)+lm1;ym]');
    drawnow
  end

endfunction

k1 = 5; k2 = 10; m1 = 0.1; m2 = 0.1; L1 = 6; L2 = 6
A = [(k1+k2)/m1 -k2/m1; -k2/m2 k2/m2];
//
function ft = f(t)
    Lambda = spec(A)
    lambda1 = Lambda(2)
    lambda2 = Lambda(1)
    omega1 = sqrt(lambda1)
    //ft = sin(omega1*t); // a)
    omega2 = sqrt(lambda2)
    //ft = sin(omega2*t); // b)
    T1 = 2*%pi / omega1;
    T2 = 2*%pi / omega2;
    ft = 0 + (sin(2*%pi*t/T1)>0); // c)
    //ft = 0 + (sin(2*%pi*t/T2)>0); // d)
endfunction
function ydot = F(t, y)
    //y=[u; v]
    B = [zeros(2,2) eye(2,2); -A zeros(2,2)];
    ydot = B*y + f(t)*[0; 0; 0; 1/m2];
endfunction
U0 = [0; 0; 0; 0]; // CI
h  = 0.05; T=10;
t0 = 0;
td = t0 : h : t0+T;

Ud = ode(U0, t0, td, F);
u1 = Ud(1,:);
u2 = Ud(2,:);
clf;
plot(td, u1, '-g', td, u2, '-r');

anime_os(td, Ud(1:2,:), m1, m2, L1, L2);