function sol=pece(system, init, period, dt)
%       [[  pece(system, init, period, 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
     % 予測 & 評価 [Adams-Bashforth method]
    x=xi;
    xdot=xdoti;
    pe=[eval(xds) eval(x2ds)];
    x=x+pe(1)*dt;
    xdot=xdot+pe(2)*dt;
    
    % 修正 & 評価 [Adams-Moulton method]
    
    ce1=[eval(xds) eval(x2ds)];
    x=xi+(pe(1)+ce1(1))*dt/2;
    xdot=xdoti+(pe(2)+ce1(2))*dt/2;
   
    % 再修正 & 評価 [Adams-Moulton method]
 %   if abs(pe(1)-ce1(1)) > 0.001
 %       ce2=[eval(xds) eval(x2ds)];
 %       x=xi+(ce1(1)+ce2(1))*dt/2;
 %       xdot=xdoti+(ce1(2)+ce2(2))*dt/2;
 %   end
    
    xi=x;
    xdoti=xdot;
    sol=[sol, [xi; xdoti]];
    t=t+dt;
    
end

