clear;

// *** 共通部分 ***
// 入力パラメータ
k = input("k = ");
d = input("d = ");
f = input("f = ");
// 初期値
x0 = input("x0 =");
dx0 = input("v0 =");
X0 = [x0; dx0];

w = %pi / 2;
t = 80;			// 時間の最大値
n = t / 0.05 + 1;	// 時間の刻み数
T = linspace(0,t,n);

// *** 常微分方程式ソルバによる数値解 ***
// 微分方程式の定義
function dx = dump(t,x)
  dx(1) = x(2);
  dx(2) = - d * x(2) - k * x(1) + f * cos(w * t);
endfunction
// 常微分方程式ソルバ
X = ode(X0,0,T,dump);

// *** グラフのプロット ***
// 時間-位置プロット
xsetech([0,0,1,0.45]);
plot(T,X(1,:),'-r');	// 数値解
xlabel("t");
ylabel("x(t)");
// 位置-速度(位相図)プロット
xsetech([0,0.5,1,0.45]);
plot(X(1,:),X(2,:),'-r');
xlabel("x");
ylabel("v");