Using RK4 method solve Dynamic 2rd PDES

Hi, I use ode45 to solve my dynamic problem well, but somehow my prof. wants me to generate the code by using RK4 method solve the problem again.I google the equations for RK4,however the result seems like going crazy. A part of code listing below：

F = @(v,u,t) M_add\(M_gr\(ForceF_rk4(t)- K_gr*u-C*v));%dv/dt G = @(v) (v);%du/dt umax = 1.5; % RK method approach; for i6 = 1:length(T)-1; k1 = h*F(v(:,i6),u(:,i6),T(i6));%dv/dt l1 = h*G(v(:,i6));%du/dt % disp(num2str(k1)); % disp(num2str(l1));

    k2 = h*F(v(:,i6)+h/4*k1,u(:,i6)+h/4*l1,T(i6)+h/4);
    l2 = h*G(v(:,i6)+h/4*k1);
%     disp(num2str(k2));
%     disp(num2str(l2));    
    k3 = h*F(v(:,i6)+1/8*h*k1+h*1/8*k2,u(:,i6)+1/8*h*l1+h*1/8*l2,T(i6)+h*1/4);
    l3 = h*G(v(:,i6)+1/8*h*k1+h*1/8*k2);
%     disp(num2str(k3));
%     disp(num2str(l3));
    k4 = h*F(v(:,i6)-1/2*h*k2+h*k3,u(:,i6)-1/2*h*l2+h*l3,T(i6)+h/2);
    l4 = h*G(v(:,i6)-1/2*h*k2+h*l3);
%     disp(num2str(k4));
%     disp(num2str(l4));
    k5 = h*F(v(:,i6)+3/16*h*k1+9/16*h*k4,u(:,i6)+3/16*h*l1+9/16*h*l4,T(i6)+h*3/4);
    l5 = h*G(v(:,i6)+3/16*h*k1+9/16*h*k4);

    k6 = F(v(:,i6)-3/7*h*k1+2/7*h*k2+12/7*h*k3-12/7*h*k4+8/7*h*k5,u(:,i6)-3/7*h*l1+2/7*h*l2+12/7*h*l3-12/7*h*l4+8/7*h*l5,T(i6)+h);
    l6 = G(v(:,i6)-3/7*h*k1+2/7*h*k2+12/7*h*k3-12/7*h*k4+8/7*h*k5);
    v(:,i6+1) = v(:,i6)+(7*k1+32*k3+12*k4+32*k5+7*k6)/90;
    u(:,i6+1) = u(:,i6)+(7*l1+32*l3+12*l4+32*l5+7*l6)/90;
end

I'm using the 2-D bar element with 10 nodes with the boundary condition of the first node did't move at all. thus the reduced mass matrix is 18*18 matrix, and the same size for stiffness matrix and damping matrix, the force is concentrated only on one node and is time depended.

the result is not even close to the result with ODE45. it going crazy to give the displacement for xx*e20 at the very beginning step.

please some one tell me what's wrong with it, and how to improve it.