clear;

// *** 解析解とソルバ解の共通部分 ***
g = 9.8;			// 重力加速度
l = g / (2 * %pi) ^ 2;		// 糸の長さ

// *** 常微分方程式ソルバによる解 ***
// 解くべき微分方程式の定義
function dx = pend(t,x)
  dx(1) = x(2);
  dx(2) = - g / l * sin(x(1));
endfunction
T = linspace(0,2,200);		// 時間ベクトル

X0 = [0; 3.1];			// 初期条件
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(TH(1,:),TH(2,:),'-r');	// プロット

X0 = [0; 6.3];			// 初期条件
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(TH(1,:),TH(2,:),'-g');	// プロット

X0 = [0; 9.4];			// 初期条件
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(TH(1,:),TH(2,:),'-b');	// プロット

X0 = [0; 12.6];		// 初期条件
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(TH(1,:),TH(2,:),'-m');	// プロット

X0 = [0; 15.7];		// 初期条件
TH = ode(X0, 0, T, pend);	// 常微分方程式ソルバ
plot(TH(1,:),TH(2,:),'-c');	// プロット

// *** グラフの体裁 ***
legend("q0 = 3.1","q0 = 6.3","q0 = 9.4","q0 = 12.6","q0 = 15.7",4);
xlabel("$\theta \mathrm{[rad]}$");
ylabel("q [rad/s]");
zoom_rect([-2,-10,8,20]);
xgrid(color(128,128,128));