２階微分方程式をルンゲクッタ法でもとめる

ode45などの関数を使用せずにルンゲクッタ法を自分で記述し２階微分方程式を解くプログラムを制作したいです。

d^2x/dt^2=-2*dx/dt-10x

という減衰振動の式を

x=x1

dx/dt=x2

とおいて刻み幅0.2でt=0~3の範囲での値を求めるプログラムです。

d^2x1/dt^2=-2x2-10x1

dx/dt=x2

という連立微分方程式にして

format long
x1=zeros(1,16);
x2=zeros(1,16);
k=zeros(1,4);
h=0.2;
kkk=[0,1/2,1/2,1,0];
a=1; x1(1,1)=1;
while a <= 15
    b=1; kk=kkk;
   while b <= 4
       k(1,b)=-2*(x2(1,a)+kk(1,b));
       kk(1,b+1)=kk(1,b+1)*k(1,b);
       b=b+1;
   end
   k(1,2)=2*k(1,2);  k(1,3)=2*k(1,3); 
   ksum=k;
   x2(1,a+1)=x2(1,a)+(1/6)*sum(ksum);
       b=1; kk=kkk;
    while b <= 4
        k(1,b)=-2*x2(1,a)-10*(x1(1,a)+kk(1,b));
        kk(1,b+1)=kk(1,b+1)*k(1,b);
        b=b+1;
    end
   k(1,2)=2*k(1,2);  k(1,3)=2*k(1,3); 
   ksum=k;
   x1(1,a+1)=x1(1,a)+(1/6)*sum(ksum);
   a=a+1;
end
t=0:0.2:3;
x=[t;x1;x2];

以上のプログラムをつくりましたがx2が速度のはずなのにずっと0で、x1も収束しません。

ちなみに完全解は

x=(1/3)*exp(-t(1,a))*(sin(3*t(1,a))+3*cos(3*t(1,a)))

[Atsushi Ueno]

format shortG

とすれば演算した値が表示されます。一つの行列xの中に絶対値の大きい値と小さい値を混在させた為、有効数字を同じ桁に表示する事が出来ていなかったのです。