matlab用解析法求二自由度阻尼系统的自由振动（先求系统状态方程

%spring2.m-用解析法求二自由度阻尼系统的自由振动（先求系统状态方程组）
m1=1;m2=9;k1=4;k2=2; 
c1=input('c1=');c2=input('c2=');  %输入阻尼系数 
x0=[1;0];xd0=[0;-1];  %初始条件
tf=50;dt=0.1;  %时间向量
M=[m1,0;0,m2];K=[k1+k2,-k2;-k2,k2];  %二阶参数矩阵
C=[c1+c2,-c2;-c2,c2];
A=[zeros(2,2),eye(2);-M\K,-M\C];  %四阶参数矩阵
y0=[x0;xd0];  %四元变量的初始条件
for  i=1:round(tf/dt)+1   %设定计算点，作循环计算
    t(i)=dt*(i-1);
    y(:,i)=expm(A*t(i))*y0;   %循环计算矩阵指数
end
subplot(2,1,1),plot(t,y(1,:)),grid
subplot(2,1,2),plot(t,y(2,:)),grid