How to solve the coupled second order differential equations using the Rung-Kutta method?

Hello,

I try to solve the two degree of freedom system(2DOF, spring-damper) with RK4.

I solve the 1DOF problem, but the answer of the above problem was not correctly.

I tried two method,

The first..

% Name: RK4_1
% Detail: Runge-Kutta 4th order examples
clear all;
close all;
h=0.001;             % [sec]
t = 0 : h: 3*pi;     % [sec]
x10(1) = 0; % [m]
x11(1) = 0; % [m/s] 
x20(1) = 0; % [m]
x21(1) = 0; % [m/s] 

 x = zeros(2, length(t));
 x(:,1) = [x10(1); x20(1)];
 v = zeros(2, length(t));
 v(:,1) = [x11(1); x21(1)];

 m1 = 1;   % [kg]
 m2 = 1;   % [kg]

 k1 = 10;  % [N/m]
 k2 = 10;  % [N/m]

 c1 = 1;   % [Ns/m]
 c2 = 1;   % [Ns/m]

 p = 10 * ones(1,length(t));   % [N]
 F = zeros(2, length(t));
 F(2,:) = p/m2;

 K11 = -(k1+k2)/m1;
K12 = k2/m1;
K21 = -k2/m2;
K22 = k2/m2;
Kmat = [K11 K12; K21 K22];

 C11 = -c1/m1;
C12 = c2/m1;
C21 = -c2/m2;
C22 = c2/m2;
Cmat = [C11 C12; C21 C22];

 f = @(t, x, y) y;
g = @(t, x, y, p) (Kmat*x + Cmat*y + p);

 for i = 1: length(t)-1
    k1 = h*f(t(i), x(:,i), v(:,i));
    l1 = h*g(t(i), x(:,i), v(:,i), F(:,i));

    k2  = h*f(t(i), x(:,i) + (1/2)*k1, v(:,i) + (1/2)*l1);
    l2  = h*g(t(i), x(:,i) + (1/2)*k1, v(:,i) + (1/2)*l1, F(:,i));

    k3  = h*f(t(i), x(:,i) + (1/2)*k2, v(:,i) + (1/2)*l2);
    l3  = h*g(t(i), x(:,i) + (1/2)*k2, v(:,i) + (1/2)*l2, F(:,i));

    k4  = h*f(t(i), x(:,i) + k3, v(:,i) + l3);
    l4  = h*g(t(i), x(:,i) + k3, v(:,i) + l3, F(:,i));

    dx = (1/6) * (k1 + 2*k2 + 2*k3 + k4);
    dv = (1/6) * (l1 + 2*l2 + 2*l3 + l4);

    x(:,i+1) = x(:,i) + dx; 
    v(:,i+1) = v(:,1) + dv;
end

 For = F';
Pos = x';
Vel = v';

 data = [Pos Vel];
xlswrite('data.xlsb',data);

And, the other method..

% Name: RK4_1
% Detail: Runge-Kutta 4th order examples

clear all;
close all;

h=0.001;             % [sec]
t = 0 : h: 10*pi;     % [sec]
x10(1) = 0; % [m]
x11(1) = 0; % [m/s] 
x20(1) = 0; % [m]
x21(1) = 0; % [m/s] 
p = 10 * ones(1,length(t));   % [N]
m1 = 1;   % [kg]
m2 = 1;   % [kg]
k1 = 10;  % [N/m]
k2 = 10;  % [N/m]
c1 = 1;   % [Ns/m]
c2 = 1;   % [Ns/m]

f1 = @(t, x1, x1d) x1d;
g1 = @(t, x1, x1d, x2, x2d) (-(k1+k2)*x1 - c1*x1d + k2*x2 + c2*x2d ) / m1;
f2 = @(t, x2, x2d) x2d;
g2 = @(t, x1, x1d, x2, x2d, F) (-k2*x2 - c2*x2d + k2*x1 + c2*x1d + F ) / m2;

for i = 1: length(t)-1
  k11 = h*f1(t(i), x10(i), x11(i));
  l11 = h*g1(t(i), x10(i), x11(i), x20(i), x21(i));
  k21 = h*f2(t(i), x20(i), x21(i));
  l21 = h*g2(t(i), x10(i), x11(i), x20(i), x21(i), p(i));

    k12 = h*f1(t(i)+0.5*h, x10(i)+0.5*k11, x11(i)+0.5*l11);
    l12 = h*g1(t(i)+0.5*h, x10(i)+0.5*k11, x11(i)+0.5*l11, x20(i)+0.5*k21, x21(i)+0.5*l21);
    k22 = h*f2(t(i)+0.5*h, x20(i)+0.5*k21, x21(i)+0.5*l21);
    l22 = h*g2(t(i)+0.5*h, x10(i)+0.5*k11, x11(i)+0.5*l11, x20(i)+0.5*k21, x21(i)+0.5*l21, p(i));

    k13 = h*f1(t(i)+0.5*h, x10(i)+0.5*k12, x11(i)+0.5*l12);
    l13 = h*g1(t(i)+0.5*h, x10(i)+0.5*k12, x11(i)+0.5*l12, x20(i)+0.5*k22, x21(i)+0.5*l22);
    k23 = h*f2(t(i)+0.5*h, x20(i)+0.5*k22, x21(i)+0.5*l22);
    l23 = h*g2(t(i)+0.5*h, x10(i)+0.5*k12, x11(i)+0.5*l12, x20(i)+0.5*k22, x21(i)+0.5*l22, p(i));

    k14 = h*f1(t(i)+h, x10(i)+k13, x11(i)+l13);
    l14 = h*g1(t(i)+h, x10(i)+k13, x11(i)+l13, x20(i)+k23, x21(i)+l23);
    k24 = h*f2(t(i)+h, x20(i)+k23, x21(i)+l23);
    l24 = h*g2(t(i)+h, x10(i)+k13, x11(i)+l13, x20(i)+k23, x21(i)+l23, p(i));

    x10(i+1) = x10(i) + 1/6*(k11+2*k12+2*k13+k14);
    x11(i+1) = x11(i) + 1/6*(l11+2*l12+2*l13+l14);        
    x20(i+1) = x20(i) + 1/6*(k21+2*k22+2*k23+k24);
    x21(i+1) = x21(i) + 1/6*(l21+2*l22+2*l23+l24);      
end

data = [x10' x11' x20' x21'];
xlswrite('data.xlsb',data);

But, both are not working correctly..

Can I get an advice for this problem?