ルンゲ・クッタ法による微分方程式の解法

以下は2階斉次線形微分方程式

d^2x/dt^2-dx/dt-2x-6=0

を

dx_1/dt=x_2
dx_2/dt=x_2+2x_1+6

のような連立1次微分方程式として解く場合の例です．

function [t,x]=runge(h)
n=0.2/h
t(1)=0
x(1,:)=[0 0]
for i=1:n
k1=h*f(t(i),x(i,:))
k2=h*f(t(i)+h/2,x(i,:)+k1/2)
k3=h*f(t(i)+h/2,x(i,:)+k2/2)
k4=h*f(t(i)+h,x(i,:)+k3)
x(i+1,:)=x(i,:)+(k1+2*k2+2*k3+k4)/6
t(i+1)=t(i)+h
end
endfunction

function dif=f(t,x)
dif(1,1)=x(2)
dif(1,2)=x(2)+2*x(1)+6
endfunction