clear;

// *** 解析解とソルバ解の共通部分 ***
// 初期振幅の入力
g = 9.8;			// 重力加速度
l = g / (2 * %pi) ^ 2;		// 糸の長さ
th0 = input("th0 = ");
k = sin(th0/2);
w = sqrt(g/l);
X0 = [0; 2 * w * k];
T = linspace(0,2,200);		// 時間ベクトル

// *** 常微分方程式ソルバによる解 ***
// 解くべき微分方程式の定義
function dx = pend(t,x)
  dx(1) = x(2);
  dx(2) = - g / l * sin(x(1));
endfunction
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(T,TH(1,:),'-r');		// プロット

// *** 解析解 ***
Tanaly = linspace(0,2,40);
THanaly = 2 .* asin(k .* %sn(w .* Tanaly, k ^ 2));
plot(Tanaly,THanaly,'og');

// *** グラフの体裁 ***
legend("o.d.e","Jacobi elliptic function",1);
xlabel("t[s]");
ylabel("$\theta \mathrm{[rad]}$");