%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% マシューの方程式
% d^2x/dt^2 + 2γ dx/dt + ω0^2 (1 + ε cos(ω t)) x = 0
% の定義
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dxdt = Mathieu(t, x, gamma, omega0, epsilon, omega)
%%%%%% 引数 %%%%%%
% t 時間
% x 状態量
%    x(1)がx(振動変位)
%    x(2)がdx/dt(振動速度)を表している
%    を表している
% gamma γ，2*γ = b/m(制動係数)
% omega0 ω0 = sqrt(k/m)(固有角振動数)
% epsilon ε
% omega ω(係数の角振動数)
%%%%%% 戻り値 %%%%%%
% dxdt xの1階微分値
%    dxdt(1)がdx/dt(= x(2))
%    dxdt(2)がd^2x/dt^2(= - 2γ dx/dt - ω0^2 (1 + ε cos(ω t)) x)
%    を表している



dxdt = [x(2); -2*gamma*x(2) - omega0^2*(1 + epsilon*cos(omega*t))*x(1)] ;
