% Runge-Kutta实现多自由度系统响应的MATLAB程序
function z=vbr_rk(rkf,u,t,x0)
%vbr_rk   vbr_rk(rkf,u,t,x0)
%       Runge-Kutta fourth order solution to a first order DE
%       t is a row vector from the initial time to the final time
%       step h.  
%       x0 is the initial value of the function.
%       The force matrix u is ordered such that the nth column of u is the force vector u evaluated at time n*dt.
%       rkf is a variable containing the name of the function file.  
%       The equations in the function file must be written in first order state space form.
%       See example vbr_rk_ex.m.
%       Example
       t=0:.005:200;                             % Creates time vector
       u=[zeros(1,length(t));sin(t*1.1)];     % Creates force matrix
       x0=[1 ;0];                             % Creates initial state vector.
%       x=vtb1_3('vbr_rk_ex',u,t,x0);         % Runs analysis.
%       plot(t,x(1,:));                       % Plots displacement versus time.
%       plot(t,x(2,:));                       % Plots velocity versus time.

n=length(t);
z=zeros(length(x0),length(t));
z(:,1)=x0;
h=t(2)-t(1);

for l1=1:(n-1);
   z1=z(:,l1);
   u1=u(:,l1);
   u2=u(:,l1+1);

   k1=h*feval('vbr_rk_ex',z1,u1,t(l1));
   k2=h*feval('vbr_rk_ex',z1+.5*k1,u1,t(l1)+.5*h);
   k3=h*feval('vbr_rk_ex',z1+.5*k2,u1,t(l1)+.5*h);
   k4=h*feval('vbr_rk_ex',z1+k3,u1,t(l1)+h);
   z(:,l1+1)=z(:,l1)+1/6*(k1+2*k2+2*k3+k4);
end

subplot(3,1,1) 
plot(t(180/0.005:200/0.005),z(1,180/0.005:200/0.005));                       % Plots displacement versus time.

subplot(3,1,2) 
plot(t(180/0.005:200/0.005),z(2,180/0.005:200/0.005));                       % Plots velocity versus time.

subplot(3,1,3) 
plot(z(1,180/0.005:200/0.005),z(2,180/0.005:200/0.005));              

%  其中 vbr_rk_ex 为动力微分方程，形如

function zd=vbr_rk_ex(zzz,u,t)
%    function for    2
%                  dx      dx
%                m --  + c -- +k x = f(t)
%                    2     dt
%                  dt                     dx
%    where m=1,k=1,c=.1, and z(1)=x, z(2)=--
%                                         dt
zd=[zzz(2);
    -.1*zzz(2)-zzz(1)+u(2)];

如有用ｒｋ方法求解微分方程组的可以参考一下．