function dx = pendubot_eq(t,x)
global a1 b1 alp3 alp2 alp5 mu2 theta1_e F

th1   = x(1);     % theta1
th2   = x(2);     % theta2
dth1  = x(3);     % dtheta1/dt
dth2  = x(4);     % dtheta2/dt

% 制御入力の計算
v = -F*(x - [theta1_e; 0; 0; 0]);

% 川田先生の解説の(22)(23)式より
ddth1 = -a1*dth1 + b1*v;
ddth2 = (alp3*dth1^2*sin(th1-th2) + alp5*sin(th2) + mu2*(dth1-dth2) - alp3*cos(th1-th2)*ddth1)/alp2;

dx = [ dth1; dth2; ddth1; ddth2 ];
      