二阶微分方程的matlab解法，以动力学方程为例

x_chuzh1=[0;0.2;0;0];
c1=0.1;
c2=0.4;
k1=2;
k2=3;
m1=4;
m2=6;
[t,x]=ode45('dyna',[0,50],x_chuzh1,[],m1,m2,c1,c2,k1,k2);
plot(t,x(:,4),'--b');
hold on
set(gca,'FontSize',10,'Fontname','Times New Roman');
function xp=dyna(t,x,flg,m1,m2,c1,c2,k1,k2)
xp=zeros(4,1);
xp(1)=x(2);
xp(2)=-(c1/m1)*(x(2)-x(4))-(k1/m1)*(x(1)-x(3));
xp(3)=x(4);
xp(4)=-(c1/m2)*(x(4)-x(2))-(k1/m2)*(x(3)-x(1))-(c2/m2)*x(4)-(k2/m2)*(x(3)-x(1))-(k2/m2)*x(3);
End