function sol=rungekutta(system, init, period, dt)
%       [[  rungekutta(system, init, dt) ]]
% この関数は二階の微分方程式をルンゲ・クッタ法により解く．
% 微分方程式は第１引数 "system"で与えられる．
% 第２引数 'period' は解析時間である．
% 第３引数 'dt' は計算の時間刻みである．
% 出力は微分方程式の変数xおよびその時間微分値xdotのデータ配列が，
% [x xdot]のフォームで出力される．
% 
% Copyright (C) 2007 清水 信行, 曽我部 潔, 小金澤 鋼一 All Rights Reserved.

% =================================================================================
% プログラムの使用条件
% 1. 本プログラムおよびデータを商用で，頒布してはならない．
% 2. 本プログラムおよびデータを著作者の許可なく営利目的で利用してはならない．
% 3. 本プログラムおよびデータの改変は，営利目的でない限り自由である．
% 4. 本プログラムおよびデータを使用することにより生じた損害などについては，著者，
%    プロクラム著作者，(社)日本機械学会，および(株)コロナ社は一切の責任を負わない．
% 5. 本プログラムの著作権は，本書編著者ならびにプログラム著作者に帰す．
% =================================================================================

xds=char(system(1));
x2ds=char(system(2));

xi=init(1);
xdoti=init(2);
t=0.;
sol=[];

while  t < period
                                        % stage 1
    x=xi;
    xdot=xdoti;
    k1=[eval(xds) eval(x2ds)];
    t=t+dt/2;
                                        % stage 2
    x=xi+k1(1)*dt/2;
    xdot=xdoti+k1(2)*dt/2;
     k2=[eval(xds) eval(x2ds)];
                                        % stage 3 
    x=xi+k2(1)*dt/2;
    xdot=xdoti+k2(2)*dt/2;
    k3=[eval(xds) eval(x2ds)];
     t=t+dt/2;
                                        % stage 4 
     x=xi+k3(1)*dt;
    xdot=xdoti+k3(2)*dt;
    k4=[eval(xds) eval(x2ds)];
     
    xi=xi+(k1(1)+2*k2(1)+2*k3(1)+k4(1))*dt/6;
    xdoti=xdoti+(k1(2)+2*k2(2)+2*k3(2)+k4(2))*dt/6;
    sol=[sol, [xi; xdoti]];
    
end
plot(sol(1,:), sol(2,:));
xlabel('x');
ylabel('xdot');
