t0=0;x0=zeros(6,1);h=0.02;
m=1*eye(3,3);   k=1*[3 -1 0;-1 2 -1;0 -1 3];

a=[zeros(3,3),eye(3,3);-inv(m)*k,zeros(3,3)];

[t,y]=ode45('m_chap2_ex2_1_sub',0,50,zeros(6,1),1e-6);

plot(t,y(:,1))