用矩阵形式表示方程：M*X''+C*X'+K*X=F(t)

方程实际上是六自由度系统的振动响应方程，就是说X是一个6×1的矩阵，M，C，K为质量，阻尼，刚度矩阵。而且C，K不是对角阵，各个自由度之间有耦合。下面是一些参量代码：

function dy=odefun(t,y)

%六自由度常微分非齐次方程

m=87.7;

I_x=2.57;I_y=4.94;I_z=3.04;

a=0.2;b=0.15;h=0.3;h_1=h/2;h_2=h/2;

k_x=1.5e5;k_y=8.75e5;k_z=1e5;

k_ox=5e4;k_oy=5e4;k_oz=5e4;

c_x=0.12*(2*sqrt(m*k_x));c_y=0.12*(2*sqrt(m*k_y));c_z=0.12*(2*sqrt(m*k_z));

c_ox=0.12*(2*sqrt(m*k_ox));c_oy=0.12*(2*sqrt(m*k_oy));c_oz=0.12*(2*sqrt(m*k_oz));

%系统参数

c_11=4*c_z+2*c_oz;

c_22=(4*c_x+2*c_ox)*b^2+(4*c_y+2*c_oy)*a^2;

c_33=4*c_x+2*c_ox;

c_44=4*c_x*h_1^2+2*c_oy*h_2^2+(4*c_z+2*c_oz)*b^2;

c_55=4*c_y+2*c_oy;

c_66=4*c_y*h_1^2+2*c_oy*h_2^2+(4*c_z+2*c_oz)*b^2;

c_16=2*c_oz*b;c_61=2*c_oz*b;

c_23=2*c_ox*b;c_32=2*c_ox*b;

c_24=-2*c_ox*b*h_2;c_42=-2*c_ox*b*h_2;

c_34=4*c_x*h_1-2*c_ox*h_2;c_43=4*c_x*h_1-2*c_ox*h_2;

c_56=4*k_y*h_1-2*k_oy*h_2;c_65=4*k_y*h_1-2*k_oy*h_2;

%阻尼矩阵元素；

k_11=4*k_z+2*k_oz;

k_22=(4*k_x+2*k_ox)*b^2+(4*k_y+2*k_oy)*a^2;

k_33=4*k_x+2*k_ox;

k_44=4*k_x*h_1^2+2*k_ox*h_2^2+(4*k_z+2*k_oz)*a^2;

k_55=4*k_y+2*k_oy;

k_66=4*k_y*h_1^2+2*k_oy*h_2^2+(4*k_z+2*k_oz)*b^2;

k_16=2*k_oz*b;k_61=2*k_oz*b;

k_23=2*k_ox*b;k_32=2*k_ox*b;

k_24=-2*k_ox*b*h_2;k_42=-2*k_ox*b*h_2^2;

k_34=4*k_x*h_1-2*k_ox*h_2;k_43=4*k_x*h_1-2*k_ox*h_2;

k_56=4*k_y*h_1-2*k_oy*h_2;k_65=4*k_y*h_1-2*k_oy*h_2;

%刚度矩阵元素

M=diag([m I_z m I_y m I_x]);

C=[c_11 0 0 0 0 c_16;

0 c_22 c_23 c_24 0 0;

0 c_32 c_33 c_34 0 0;

0 c_42 c_43 c_44 0 0;

0 0 0 0 c_55 c_56;

c_61 0 0 0 c_65 c_66];

K=[k_11 0 0 0 0 k_16;

0 k_22 k_23 k_24 0 0;

0 k_32 k_33 k_34 0 0;

0 k_42 k_43 k_44 0 0;

0 0 0 0 k_55 k_56;

k_61 0 0 0 k_65 k_66];

%质量矩阵，刚度矩阵，阻尼矩阵；

F=[1*sin(10*t);2*sin(10*t);3*sin(10*t);4*sin(10*t);5*sin(10*t);6*sin(10*t)];

%激振力列阵

x=y(1:6);%前六项对应6个y，也就是位移

dx=y(7:12);%后六项对应6个y的一阶导数，也就是速度

dy=[x;M\((F-K*x-C*dx))];

tspan=[0 0.02];%微分区间

y0=zeros(12,1);%初值共12项，即初始位移和速度，我这里默认全部为0

[t,y]=ode45(@odefun6,tspan,y0);

plot(t,y)