clear;

// *** 計算条件 ***
g = 9.81;	// 重力加速度
k = 0.8;	// 反発係数
// 時間
ts = 0;		// 開始時刻
te = 10;	// 終了時刻
// 初期条件
x0 = 10;	// 初期位置
v0 = 15;	// 初速度

// *** 常微分方程式の定義 ***
function dx = fall(t,x)
  dx(1) = x(2);	// dx/dt = v
  dx(2) = -g;	// dv/dt = -g
endfunction

// *** 解くべき非線形方程式 ***
function x = bound(t)
  X = ode([x0; v0], ts, t, fall);
  x = X(1);
endfunction

// *** メイン ***
// 初期化
tb = ts;	// 衝突時刻
X = [];		// プロット用の位置と速度
T = [];		// プロット用の時間
// 終了時刻まで繰り返し
while tb < te
  // 衝突時刻を計算
  tb = fsolve(2 * te, bound);
  // プロット用の時間ベクトルを作成
  if tb > te then
    // 終了時刻まで
    time = linspace(ts, te);
  else
    // 衝突時刻まで
    time = linspace(ts, tb);
  end
  T = [T, time];
  // プロット用の位置と速度を計算
  X = [X, ode([x0; v0], ts, time, fall)];
  // 次の繰返しのための初期条件
  x0 = 0;		// 初期位置
  v0 = -k * X(2,$);	// 初速度
  ts = tb;		// 計算開始時刻
end

// *** グラフのプロット ***
// 位置のプロット
subplot(2,1,1);
plot(T, X(1,:),'g');
xgrid(color("gray"));
xlabel("Time");
ylabel("Position");
// 速度のプロット
subplot(2,1,2);
plot(T, X(2,:),'r');
xgrid(color("gray"));
xlabel("Time");
ylabel("Velocity");