clear;

// ************************************
// 
//               共通部分
// 
// ************************************
// *** 入力パラメータ ***
a = 1;
b = input("b = ");
k = input("k = ");
// 時間ベクトル
T = linspace(0,10);

// ************************************
// 
//                数値解
// 
// ************************************
// 解くべき微分方程式
function dx = balanced(t,x)
  dx = k * (a - x) * (b - x);
endfunction
// 初期値
x0 = 0;
// 微分方程式の数値解
X = ode(x0,0,T,balanced);

// ************************************
// 
//                解析解
// 
// ************************************
g = (a - b) / 2 * k;
if a == b then
  Xa = a - 1 ./ (k .* T + 1 / a);
else
  Xa = a * b * (exp(g * T) - exp(- g * T)) ./ (a * exp(g * T) - b * exp(- g * T));
end

// ************************************
// 
//               プロット
// 
// ************************************
plot(T,X,'or');
plot(T,Xa,'-g');
xtitle(strcat(["Chemical Reaction: ","a = ",string(a),", b = ",string(b),", k = ",string(k)]));
xlabel("$t [s]$");
ylabel("$N_c$");
legend(["Numerical","Analytical"],4);